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THE FIRST GALAXIES: ASSEMBLY UNDER RADIATIVE FEEDBACK FROM THE FIRST STARS 
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ABSTRACT 

We investigate how radiative feedback from the first stars affects the assembly of the first dwarf 
galaxies. To this end we perform cosmological zoomed smoothed particle hydrodynamics simulations 
of a dwarf galaxy assembling inside a halo reaching a virial mass ~ 10 9 M Q at z = 10. The simulations 
follow the non-equilibrium chemistry and cooling of primordial gas and the subsequent conversion of 
the cool dense gas into massive metal-free stars. To quantify the radiative feedback, we compare a 
simulation in which stars emit both molecular hydrogen dissociating and hydrogen/helium ionizing 
radiation with a simulation in which stars emit only molecular hydrogen dissociating radiation, and 
further with a simulation in which stars remain dark. Photodissociation and photoionization exert 
a strong negative feedback on the assembly of the galaxy inside the main minihalo progenitor. Gas 
condensation is strongly impeded, and star formation is strongly suppressed in comparison with the 
simulation in which stars remain dark. The feedback on the gas from either dissociating or ionizing 
radiation implies a suppression of the central dark matter densities in the minihalo progenitor by 
factors of up to a few, which is a significant deviation from the singular isothermal density profile 
characterizing the dark matter distribution inside the virial radius in the absence of radiative feedback. 
The evolution of gas densities, star formation rates, and the distribution of dark matter becomes 
insensitive to the inclusion of dissociating radiation in the late stages of the minihalo assembly, and it 
becomes insensitive to the inclusion of ionizing radiation once the minihalo turns into an atomically 
cooling galaxy. The formation of a rotationally supported extended disk inside the dwarf galaxy is 
a robust outcome of our simulations not affected by the inclusion of radiation. Our estimates of the 
observability of the first galaxies show that dwarf galaxies such as simulated here will be among the 
faintest galaxies the upcoming James Webb Space Telescope will detect. 

Subject headings: cosmology: theory - galaxies: formation - galaxies: high-redshift - stars: formation 
- hydrodynamics - radiative transfer 



1. INTRODUCTION 

The birth of star-forming galaxies a few hundred mil- 
lion years after the Big Bang marks an important mile- 
stone in the history of our universe. The spectacular 
images returned by the Hubble Space Telescope (HST) 
in the past few years have already allowed us to probe 
into the first billion year of the universe. The last few 
years have also seen the development of an observation- 
ally testable theory of the formation of the first galaxies. 
As ongoing and upcoming observations, such as with the 
James Webb Space Telescope (JWST), are about to push 
to ever earlier times, approaching the epoch of the very 
first galaxies, this theory i s put to ever more stringent 

tests (for a review see, e.g., lDunlodl20Tl. 

Both analytica l arguments (e.g.lTegjmark et al.l 119971 : 
iNaoz et all 120061: iTseliakhovich et all 1201 ID and sim- 
ulations (e .g.. lAbel et al.l l2002t iBromm et al.l l2002t 
L312Q1K); 



iStacv et"al~ll2011[ ) suggest that the first stars have formed 
at z > 30 inside dark matter minihalos with virial 
temperatures > 10 3 K, correspondi ng to halo masses 
~ 10 5 — 10 6 M Q (for a review, see IBromm fc Larson! 
2004). The metal- free primordial gas inside these mini- 
halos cools and condenses to reach the high densities 
needed to form stars primarily through the radiative de- 
excitation of rovibra tionally excited molecular hydrogen 
(e.g.. lHaiman et al.l [l996bl lAbel etail 119971 ). As the 
minihalos grow in mass, their virial temperatures in- 
crease and, after reaching ~ 10 4 K, become sufficiently 



large to collisionally excite atomic hydrogen. This gives 
birth to the first atomically cooling galaxi es with typical 
masses ~ 10 7 - 10 s M^ at z > 15 (e.g.. lOh fc Haimanl 
200 21 IWise fc Abelll2ll07b1: IGreif et al.l l2008: for a review 
see lBromm fc Yoshidal l2011'). The first atomically cool- 
ing galaxies then evolve into the first dwarf galaxies with 
characteristic masses > 10 9 M (e.g. JMashchenko et all 
120081 IWise et al.ll2012l ). The aim of the current work is 
to investigate the assembly of such galaxies under the 
radiative feedback from the first stars. The mass-scale 
marked by the first dwarf galaxies is closely related to a 
number of key open issues, some of which are outlined 
below. 

Suppression of star formation by radiative feedback. 
The radiation emitted by the first stars has a pro- 
found effect on subsequently fo rming stars and galaxie s 
(for a comprehensive overview. iCiardi fc Ferraral 120051) . 
Radiation in the Lyman-Werner (LW) bands dissoci- 
ates molecular hydrogen, the main coolant in minihalos 
(e.g.. lHaiman e t al. 1997). Hydrogen-ionizing radiation 
heats the gas inside the first halos and the intergalac- 
tic medium (IGM). The associated increase in pressure 
drives the gas outside halos with virial temperatures < 
10 4 K, suppressing star formation in both minihalos and 
the first atomic cooling halo s (e.g., iThoul fc Weinberg] 
119961 : IBarkana fc Loebl I1999D . The increased pressure 
in the IGM impedes the accretion of gas onto these 
low-mass halos, an effect known as Jeans-filtering (e.g., 
IShapiro et all 119941: iGnedin fc Huil H998: Oka moto etaLl 
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2008). In contrast, star formation inside the first dwarf 
galaxies should be more robust to this negative radiative 
feedback as their deeper gravitational potentials allow 
them to hold on to their gas more strongly. 

The sources of reionization. The first galaxies are 
thought to have started the reionization of the universe, 
which is the transformation of the cosmic hydrogen from 
its early neutral to its present ionized state that occurred 
during the first billion year after the Big Bang (for re- 
views see, e.g., Bar kana fc Loebl 120011 : iFurlanetto et al.l 
|2(M iMeiksirJ 12009^ 7" However, whether galaxies could 
sustain reionization and drive it to completion is a 
question of significant debate (e.g.. iBouwens et aLll2012l : 
iFinkelstein et alJl2012D . The largest uncertainties are re- 
lated to our poor knowledge of the escape fraction, i.e., 
the fraction of ionizing photons that leave the galaxies 
unabsorbed and are thus available to reionize the IGM, 
and of the abundance and ionizing luminosities of low- 
mass gala xies too faint to be d etected in current sur- 
veys (e.g., iRobertson et alJfeOlOl) . The most recent de- 
terminations of the UV luminosity density at z > 7 sug- 
gest that a significant contribution from faint, yet to be 
observed, low-mass galaxies is likely needed to sustain 
reion izati on in the recombining gas ()Finkelstein et all 
[2011 see lKuhlen fc Faucher-Giguerell2012l for a compre- 
hensive discussion). Because of their increased robust- 
ness against stellar feedback, dwarf galaxies are among 
the low-mass galaxies expe cted to be especially efficient 
sources of reionizati on (e.g. lChoudhury fc Ferraral 120071 : 
iRaicevic et al.l[20TTI ). 

The origin of the Milky Way (MW) satellites. Dissi- 
pationless simulations of dark matter subhalos around 
MW-like galaxies imply the existence of a large popu- 
lation of satellite galaxies only few of which are cur- 
rently observed. A number of solutions have been offered 
to explain this missing satellite problem (iMoore et al. 
19991: IKlvpin et al.l [l999l fo r reviews see, e.g., iBullockl 



20iq iRicottil I2010t iMaverl [2010h . including the sup- 



pression of star formation in low-m ass halos by stel- 
lar feedback and reionization ( e.g.. iRicotti fc Gnedinl 



200^ iSalvadori fc Ferraral [20091 iBovill 



Sawala et al. 2012: Brooks 



^Ricottil l200£. 
Zolotovll2012h . the trans- 



formation of low-mass halos by tidal forces and ram- 
pressure strip ping upon their e ntry in the MW virial 
region (e.g., iMaver et al.l [2007), modifications of the 
cold dark matter structure formation paradigm (e.g., 
Lovel et al.l l2012fl. and others (e.g., Bow fc Dvorkinl 



20121 : IVera-Oiro et al.ll2l!H IWang et al.ll2012D . However. 



current theories still struggle to explain the abundance 
and properties of the observed MW satellites across the 
luminosity range, from the recently discove red ultra-faint 
to the long-known classical satellites (e.g. | Strigari et al.1 
[20081: IBovill fc RlcottJ[20TTI: IBovlan-Kolchin et alJl201lT r 
The most massive of the simulated MW dark matter sub- 
halos have progenitors with masses 10 8 — 10 10 M Q at 
z > 6 (e.g.. IBovlan-Kolchin et al1l2012|) . This suggests 
an intimate relation with the first dwarf galaxies, and 
renders investigations into these objects a promising tool 
to understand the origin of structure in galaxies such as 
the MW. 

The first disk galaxies. Simulations of the first atom- 
ically cooling galaxies, i.e., galaxies inside halos with 
masses 10 7 — 10 8 M at redshifts z ~ 15, r eveal a highly 
irregular morphology of the halo gas (e.g.. iWise fc Abell 



I2007bb l&iffet~aTl [20081 IRegan fc Haehneltl [20091) . On 
the other hand, simulations of galaxies inside halos with 
larger masses > 10 9 M Q at lower redshifts z < 6 of- 
ten find the halo gas organized in rotat i onally sup- 
porte d disks (e.g., iMashchenko et al.1 120081 : IWise et all 
2012). These findings suggest a dwarf-size mass scale 
~ 10 8 — 10 10 Mq for the transition to disk-like mor- 
phologies, and the emergence of the first disk galaxies 
at z > 6. Physical processes to imprint such a scale 
include the turbulence generated by the cold inflow of 
gas along filaments which characte rizes gas accretion by 
the first atomic cooling halos (e.g ., IWise fc Abell [2007b: 
IWise et all 120081: iGreif et aLl l2008lh and stellar feedback 
(e.g.. iKaufmann et al. Whether the first halos 

may host disks is an important open issue, affecting es- 
timates of, e.g., the escape of ionizing photons into th e 
IGM (e.g.. IGnedin"eraT]^00a [Conrov fc Kratterll2012ft . 
or the abili ty of massive black holes to accrete gas and 
grow (e.g.. lEisenstein fc Loebl 119951: iKoushiappas et al.1 
120041: lLodato fc Natarajanll2006HPetri et al.ll2012| ). 

The faintest galaxies JWST will see. The faintest 
z > 6 galaxies HST has so far revealed have esti- 
matcd stellar masses > 10 8 MfD (ILabbe et all 120101: 
Fin kelstein et al.ll2010l: (Curtis-Lake et al.ll2012t)~ Future 
observations with upcoming telescopes such as the JWST 
will allow to search for galaxies down to still lower stel- 
lar masses and out to higher redshifts, thus promising 
to test our theories of the formation of the first stars 
and galaxies. However, most studies agree that even 
JWST will not be sensitive enough to detect the stel- 
lar radiation emitted from inside the minihalos and the 
first atomic co ol ing halos (e.g. . lHaiman fc Loebl 119981: 
Oh et all 120011: Uohnson et al.1 120091: iZackrisson et atl 



2012tlRvdberg et al.ll2012l ). JWSTm&v detect the stellar 



radiation from so me of these objects if th ey are gravita- 
tional lensed (e.g.. IZackrisson et al.| [2012l . But most of 
the stellar light collected by JWST from high redshifts is 
expected to come from dwarf galaxies mor e massive than 
the first atomically cooli ng galaxies (e.g.. Uohnson et"aTI 
2009; lPawlik et alj|201ll) . 

Motivated primarily by the exciting prospects for ob- 
servations with the upcoming JWST, we have previously 
presented cosmological simulations of a dwarf galaxy 
assembling inside a halo reaching 10 9 M Q at z = 10 
(jPawlik et al.ll20lTI ) . The simulations were performed us- 
ing the Smoothed Particle Hydrodynamics (SPH) tech- 
nique and achieved high resolution by zooming in a select 
region around the galaxy. Following the non-equilibrium 
chemistry and cooling of primordial gas, the simulations 
tracked the evolution of the dwarf galaxy starting from 
before its birth inside a minihalo. An intriguing outcome 
was the formation of a rotationally supported extended 
disk just prior to the final simulation redshift. However, 
our previous simulations did not account for star forma- 
tion and the associated feedback. As explained above, 
stellar feedback has the potential to significantly affect 
the assembly of the gas inside low-mass halos. 

In this study we present a new set of simulations sim- 
ilar to our previous simulations, but extending them by 
including star formation and radiation. We focus on the 
radiative feedback from LW and ionizing radiation. 1 To 

1 We will use the terms LW radiation and dissociating radiation 
interchangeably. 
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judge the impact of radiative feedback we will compare 
a simulation that includes both LW and ionizing radia- 
tion with a simulation that includes only LW radiation 
and further with a simulation in which no radiation is 
emitted. Note that the simulations do not account for 
supernova (SN) feedback or metal enrichment, a limita- 
tion we will discuss in Section [7] below. The simulations 
are designed primarily to address the impact of radiative 
feedback on the assembly of the emerging dwarf galaxy, 
and on the formation of galactic disks inside it. However, 
we will also briefly discuss the impact of radiative feed- 
back on the assembly of galaxies in the neighborhood of 
the simulated dwarf galaxy. Our simulations enable us 
to provide an improved estimate of the observability of 
the first galaxies with JWST. 

The organization of this paper is as follows. In Sec- 
tion [2] (as well as in the appendix) we describe our nu- 
merical techniques, and in Section [3] we describe the set 
of simulations that we have carried out. In Sections 2] 
and [5] we present the results of our simulations, subse- 
quently discussing the assembly of the dwarf galaxy and 
the radiative feedback on the IGM and the neighboring 
galaxies. In Section [6] we use the simulated star forma- 
tion rates to estimate the observability of the first galax- 
ies with JWST. In Section [7] we discuss our results and 
also address some of the most important limitations. In 
Section [51 we summarize our work. 

Throughout this work we assume the ACDM cos- 
mological model with parameters Q m — 0.258, tth — 
0.0441, 0\ = 0.742, erg = 0.796, n s = 0.963, and h = 
0.719, which are consistent with the most recent anal- 
ysis of the observations wi th the Wilkinson Mic rowave 
Anisotropy Probe satellite (|Komatsu et al.l 120101) . Dis- 
tances are expressed in physical (i.e., not comoving) 
units, unless noted otherwise. We will make use of the 
species number density fractions with respect to hydro- 
gen rj a = n Q /riH, where a labels the chemical species. 

2. NUMERICAL METHODS 

In this section we describe the numerical techniques 
employed. The simulations presen ted below are identi- 
cal to the simulations described in iPawlik et al.l (| 2 1 If ) . 
except for the inclusion of star formation and dissoci- 
ating and ionizing radiation. We will therefore only 
briefly review the s imulation techniques already used in 
IPawlik et al.l ([201 ID . We will focus on the description of 
the techniques used to model the formation of stars and 
the dissociative and ionizing impact of stellar radiation 
on the gas. We remind the reader that the simulations 
do not account for SN feedback or metal enrichment, as 
discussed further in Section [7J 

2.1. Gravity and Hydrodynamics 

We use a modified version of the iV-body/TreePM 
Smoothed Part i cle Hydrodynamics ( SPH) code GADGET 
([Springe!! 120051: iSpringel et aL1[200ll: iSchave et al1[20Tol) 
to perform a suite of zoomed cosmological hydrodynam- 
ical simulations of the assembly of a halo that reaches 
a virial mass M V - 1T ss 10 9 M Q at redshift z — 10. The 
simulations are initialized at redshift z — 127 in a box of 
size 3.125 h^ 1 Mpc comoving. Initial particle positions 
and velocities are obtained by applying the Zeldovich ap- 
proximation (| Zeldovichl I1970T) to particles arranged on a 



Cartesian grid. We adopt a transfer function for mat- 
ter perturbations generate d with CMBFAST (version 4.1; 
ISeliak fc Zaldarriagal[l996l) . 

We first perform a low-resolution simulation without 
star formation down to redshift z = 10, and locate the 
most massive halo, with virial mass M w -„ ps 10 9 M© and 
virial radius r v - u = 3.1 kpc. We then trace the particles 
found within 3r v j r from the most bound particle of this 
halo back to their locations at the start of the simula- 
tion. We hierarchically refine the initial particle setup 
using a nested sequence of cubical patches ("zooms") 
centered on the traced particles, inside which we increase 
the mass resolution by successive factors of 8 with each 
increasing level of zoom. In the zoom with the high- 
est mass resolution, which entirely contains the traced 
particles and which we refer to as the refinement region, 
gas (dark matter) particles have masses m g = 484 M© 
(ttidm = 2350 Mq). The simulations are performed 
with the gravitational forces softened over a sphere of 
Plummer-cquivalent radius e = 0.1 h^ 1 kpc comoving 
applied to all particles. 

The particle dynamical variables, such as position, ve- 
locity, and density, are evolved in time using individ- 
ual particle gravito-hydrodynamical time steps deter- 
mined by the smaller of the dynamical time step and 
the C ourant time step (e.g., Equation 16 in iSpringell 
2005), which is the standard GADGET time stepping 
scheme. We do not explicitly limit the particle gravito- 
hydrodynamical time steps by either the chemical time 
or the radiative cooling or radiative heating time. How- 
ever, the chemistry, radiative cooling and radiative heat- 
ing described below are solved by subcycling the small- 
est gravito-hydrodynamical time step among all parti- 
cles in th e sim ulation on the relevant time scales (see 
Appendix I A. 4j) . 

2.2. Chemistry and Cooling 

We assume that the gas is of primordial composition 
with hydrogen mass fraction X = 0.75 and a helium 
mass fraction Y = 1 — X. We use a modified ver- 
sion of the implicit solver dvode ([Brown et al.1 119891 ) 
to follow the non-equilibrium chemistry and cooling of 
H 2 , D, HD, D+, H+, H, D, and He,' and we include 
H and Hj" assuming their collis i onal equilibrium abun- 
dances ([Johnson fc BromirJ|2Tjol iGreif et al.l[2010l ). We 
consider all relevant radiative cooling processes: cooling 
by collisional ionization, collisional excitation of atomic 
and molecular lines, the emission of free-free and recom- 
bination radiation, and Compton cooling by the CMB. 
Once stars form and emit radiation, the chemical and 
thermal evolution of the gas is also affected by the pho- 
todissociation of molecular hydrogen and deuterium, and 
photoionization of hydrogen and helium, as we describe 
in Sections 12.51 and 12.61 below. It should be kept in mind 
that at high gas densities, a Jeans floor employed to avoid 
artificial fragmentation artificially increases the gas tem- 
perature and affects th e distribution and d ynamics of the 
gas (see Equation 1 in IPawlik et alJl2TjTTl ). 

2.3. Star Formation 

Star formation is a complex astrophysical phe- 
nomenon, many details of w hich remain to be under - 
stood (for a review see, e.g., iMcKee fc Ostrikerl l200l . 
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In the MW and in nearby galaxies, star formation is ob- 
served to occur inside a hierarchy of clouds with masses 
~ 10 4 — 10 6 Mq. The clouds are transformed into stars at 
a rate = /O c i/r*, where = rg/ eg is the star formation 
time scale, Tg — [3n/ (^Gpd)} 1 ^ 2 the free fall time at the 
characteristic density p c \ of the star-forming clouds, and 
eg the star formation efficiency per free fall time. Star 
formation is observed to be a slow process with a typical 
efficiency per free fall time of only eg ~ 0.01, indepen- 
dent of the characteristic densities of the star-forming 
clouds in the rang e rend = PciX/mn ~ 10 — 10 5 cm~ 3 
(|Krumholz fc TarJl2007h . 

While gas masses ~ 10 4 — 10 6 M Q that characterize 
star-forming clouds in the nearby universe are resolved, 
our simulations lack the resolution and physical detail to 
follow the formation of individual stars. Hence, we can- 
not exploit our simulations to estimate the star formation 
rates (SFRs) inside individual star-forming clouds from 
first principles. Instead, we adopt a phenomenological 
model that specifies how quickly gas turns into stars. The 
model is motivated by and consistent with investigations 
of star formation in the nearby universe. We note that 
prescriptions for treating star formation in simulations of 
high-rcdshift galaxies are commonly calibrated with re- 
lations obtained from the local universe. This practice is 
necessitated by the current lack of direct observations of 
star formation at high redshifts. The phenomenological 
approach is sufficient to allow us to investigate the effect 
of stellar dissociating and photoionizing radiation on the 
interstellar gas and the IGM. 

We restrict star formation to occur only in regions 
with gas densities exceeding a threshold density, hh > 
nsp, where we set ngF = 500 cm~ 3 . This is some- 
what lower than the densities ~ 10 4 cm~ 3 of metal- 
free clouds on the v erge of collapse t o form stars (e.g., 
iBromm et all l2002t I Abel et al.ll2002f ). It is also lower 
than the similar densities of metal-free clouds able to 
shield their H2 from external dissociating radiation (e.g., 
iSafran ck-Shra der et alj [2012f) . Because our simulations 
focus on the assembly of galaxies inside more massive 
halos, we must adopt, for reasons of computational via- 
bility, a lower resolution than employed in these previous 
works, preventing us from adopting still larger star for- 
mation threshold densities. 

We adopt a star formation time scale that depends 
solely on the threshold density for star formation nsF, 



Tff(nsF) 



201.2 



eg 



Myr ( 



500 cm- 



-1/2 



(1) 



where we have set eg = 0.01. Hence, the star formation 
time scale is independent of the gas density nn > nsp. 
This amounts to assuming that all star formation occurs 
inside clouds with characteristic density nn.ci = Hsf, and 
is consistent with both our limited resolution and the ob- 
servations in the nearby universe. A similar star forma- 
tion recipe using a density-inde pendent s t ar for mation 
time has been employed in, e.g.. iKravtsovl ()2003[ ). 

Our numerical imple mentation of the star formatio n 
law is identical to that of Schave fc Dalla Vecchial {2008). 
The star formation law is interpreted stochastically, and 
the probability that a star-forming gas particle is turned 
into a star particle in a time interval At is given by 
min(At/r Jr , 1). Individual gas particles give rise to a 



single star particle, i.e., the masses of the star par- 
ticles are identical to that of the gas particles from 
which they are formed. We impose an upper limit 
= max(1.5 x 10 4 K, 10 5 x T floor ) on the tempera- 
ture at which gas is allowed to form stars, where Tg OOI is 
the minimum temperature set by the Jeans floor. 

2.4. Population Synthesis 

We interpret the star particles in our simulations 
as simple stellar populations, i.e., instantaneous stellar 
bursts that are characterized by an initial mass func- 
tion (IMF), metallicity, and age. We compute the time- 
dependent hydrogen and helium ionizing luminosities 
Qhi, QhcI, and QhcII, as well as the luminosities Qlw in 
the LW band with energies 11.2 — 13.6 eV, of these star 
form ation bursts usin g the population synthesis models 
from iSchaererl (|2003| ). The m odels ass u me a power- law 
IMF p(m*) ex m~ a with the ISalpeterl (|1955|) exponent 
a = 2.35 but allow for a variation of the range of the 
stellar masses sampled from the power- law distribution . 
We describe the stellar bursts using the ISchaererl (pOOl ) 
zero metallicity models with initial masses in the range 
50 — 500 M Q . The age of a burst is the time difference 
between the simulation time at which the star particle 
was created and the current simulation time. For ref- 
erence, the luminosities at zero age main sequence are 
Qhi = 10 47 - 98 , QhcI - 10 47 - 80 , QhcII - 10 4705 , and 



Qlw = 10 photons Mq s . We only consider star 
particles inside the refinement region, and we do not fol- 
low the propagation of radiation outside this region. 

2.5. Photodissociation 

Molecular hydrogen H2 and deuterated hydrogen HD 
are photodissociated upon absorption of radiation in the 
LW band belonging to the photon energy range 11.2 eV— 
13.6 eV. The r ate of photodissociation of H2 is (e.g., 
lAbel et al.iri997h 



fc LW = l.l x 10 8 



erg Hz 1 s 1 cm 



= 1.38 x 10~ 12 s -1 J 2 i, 



(2) 



the characteris- 
tic LW flux, 



where z^lw = 12.4 eV is 
tic LW frequency, ^(z/lw) is 

J21 = J,/(W 21 erg Hz-' 
normalized LW mean intensity, and J w = F v (yuf/) / A-k . 
We set the rate for photod issociation of HD iden- 
tical to that of H2 (e.g., Glover fc Jappsenl 120071 : 



IWolcott-Green fc Haimar3l20lC 

We write the LW intensity J21 as_the sum of the cos- 
mological LW background intensity J21 and the local LW 
intensity J^ produced by the stellar populations repre- 
sented by the star particles in the refinement region of 
our simulations, J%\ = J21 + J\°i '■ The refinement region 
is too small to follow the build-up of the cosmological 
LW background in a self-consistent manner. We there- 
fore treat the intensity of the LW background as a free 
parameter, and approximate its evolution using 



(3) 



where J2i,o is the intensity of the LW background at z . 
Setting J2L0 = 1 and zq = 10, Equation (|3|) provides 
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a good fit to the LW backg r ound evolution presented in 
Figure 1 of lGreif fe Brom m ( 2006) , and is also consistent 
with computations of the LW background in other works 
(e.g. JWise fc Abe]ll20"o"3 lAhn et al.ll2009l ). 

For reasons of computational efficiency we determine 
the contributions from the AT* star particles in the refine- 
ment region to the intensity jj° c — Ylj=i Jj evaluated 
at the location of gas particle i in the optically thin ap- 
proximation (e.g.. IWise fe A bel 2008a|), i.e., 



J* = ■ 
3 AttSvlw 47rr?- 



or, in units of 10 erg Hz s cm sr 

\ ( m i 



21,j ~ 1 



aw,j 



10 47 M^ 1 s- 1 ) \500 M Q J \ 1 kpc 



(4) 



(5) 



Here, Qlwj is the photon luminosity of star particle 
j per unit mass in the LW frequency band, (5^lw = 
(13.6 eV - 11.2 eV)/h P is the width of the LW band, h P 
Planck's constant, m* the mass of star particle j, and r%j 
the distance between the gas and the star particle. 

For simplicity, we approximate the time-dependent LW 
luminosities Qlw computed in Section \2A\ by their zero 
age main sequence values, and we assume that the stars 
emit LW radiation for 3 Myr (jSchaererl 120031 ) . We ap- 
proximately account for radiative transfer (RT) effects 
by attenuating the total LW inte nsity J21 by a local 
self-sh ielding factor (Equation 10 in lWolcott-Green et al.l 
1201 lbl with a = 1.1). The self-shielding factor de- 
pends on the column density of molecular hydrogen, 
which we compute in the local Jeans approximation (e.g., 
IShang et al]l20i0h . 

When computing the contribution to the LW inten- 
sity from local sources using Equation ((4]), we ignore 
the absorption of LW radiation along the way to the 
absorbing gas particle. Near emitters, this is justi- 
fied because molecular hydrogen is efficiently collision- 
ally d issociated inside th e H II regions surroundin g 
them (| Johnson et alj|2007t but see lRicotti et all 120021 ). 
Outside the H II regions, this approximation is good 
as long as the intergalactic molecular hydrogen frac- 
tion does not substantially exceed its primordial value, 
r /H 2 ~ 10~ 6 . In this case, the optical depth for ab- 
sorption in the LW band remains insignificant out to 
distances > 1 Mpc co moving (e.g.. iRicotti et al.l I2001L 
their Figure 11; see als o IGlover fe Bra nd 2003, their Fig- 
ure 7; lAhn et a l. 2009) much larger than the size of the 
refinement region. Finally, we also ignore the absorption 
of LW photons by atomic hydrogen series lines, which 
also i s not significant on t h e small sca l es sim ulated here 
(e.g., lHaiman et all [20001: lAhn et al.l [20091 ) . Our im- 
plementation of ph otodissociation is similar to that in 
Uohnson et al.ll2012[ but we do not account for the photo- 
detachment of H~. 

2.6. Photoionization 

In this section we describe our implementation of pho- 
toionization by stellar radiation. This implementation 
involves two main steps. First, we transport ionizing 
photons radially from the star particles through the sim- 
ulation box, and compute the fraction of the photons ab- 
sorbed by atomic hydrogen and helium. Second, we infer 



the associated photoionization and photoheating rates. 
Here we will only give a brief overview. We present a de- 
tailed description of the implementation of the RT and 
the computation of the photoionization and photoheat- 
ing rates in the Appendix. There we will also discuss 
tests of this implementation. The coupling of the RT 
with the hydrodynamical evolution is achieved by passing 
the photoionization and photoheating rates, along with 
the photodissociation rates described in the previous sec- 
tion, to the non-equilibrium solver for the chemical and 
thermal evolution of the gas described i n Sec tion l2.21 and 
this coupling is described in Appendix IA. 41 

We transport the ionizing radiation emitted by 
the star p articles using the multi- fre quency RT code 
traphic (jPawlik fe Schave I 120081: iPawlik fe Schavd 
120111) . TRAPHIC solves the time-dependent RT equation 
by tracing photon packets emitted by source particles 
at the speed of light and in a photon-conserving man- 
ner through the simulation box. The photon packets are 
transported directly on the spatially adaptive, unstruc- 
tured grid traced out by the SPH particles, which allows 
one to exploit the full dynamic range of the SPH simula- 
tions. A directed radial transport of the photon packets 
from the sources is accomplished despite the irregular dis- 
tribution of SPH particles by guiding the photon packets 
inside cones. A photon packet merging technique ren- 
ders the computational cost of the RT independent of 
the number of ionizing sources. The transport of pho- 
tons is discretized in RT time steps At T , after each of 
which the chemical and thermal evolution of the gas is 
advanced based on the number of absorbed photons. 

We make the following approximations specific to the 
current work. Ionizing photons with energies above 
13.6 eV are transported using a single frequency bin, and 
the absorption of the photons by neutral hydrogen and 
neutral and singly ionized helium is computed in the grey 
approximation. Each star particle emits photon packets 
to its neighboring SPH particles once per RT time step in 
a set of tessellating emission cones centered around 8 dif- 
ferent directions. The effective angular sampling of the 
surrounding volume is larger than implied by this num- 
ber of directions thanks to the splitting of photon packets 
among neighbors inside the same emission cone, and be- 
cause of the randomization of the emission directions at 
each RT time step. The photons are transported radially 
away from the star particles by tracing them downstream 
inside transmission cones with solid angle 47r/128. This 
angular resolution is sufficiently high to track the delay 
of the ionization fronts around individual halos by dense 
filaments, giving rise to the typical "butterfly" shape of 
ionized regions, as shown in Figure [TJ To reduce the com- 
putational cost, we limit the propagation of photons to 
at most a single inter-particle distance per RT time step, 
which approximates the full time-dependent RT in the 
limit of small RT time steps. These approximations are 
discussed in further detail in the Appendix, in particular 
in Appendix IB. 31 

3. SIMULATIONS 

We employ a set of three simulations to study the 
effects of radiative feedback on the assembly of high- 
redshift galaxies. Simulation LW+RT includes both star 
formation and the LW and ionizing radiation emitted 
by the stars, as well as an imposed LW radiation back- 
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TABLE 1 
Simulation parameters. 



Simulation m gas a m DM b e c n SF d LW C RT* 

LW+RT 4.84 X 10 2 2.35 X 10 3 0.1 500 yes yes 

LW 4.84 x 10 2 2.35 x 10 3 0.1 500 yes no 

NOFB 4.84 x 10 2 2.35 x 10 3 0.1 500 no no 



a Gas particle mass in the refinement region (Mq). 
b Dark matter particle mass in the refinement region (Mq). 
c Gravitational softening radius (h -1 kpc comoving). 
d Star formation threshold density ngp (cm -3 ). 
° Photodissociation by stellar LW photons and the LW background. 
Radiative transfer of stellar ionizing photons. 



ground, as described in Section [2] We will focus on dis- 
cussing results of this simulation. We will often present 
our results by comparing this simulation with simulation 
LW, which is identical except that the emission of ion- 
izing radiation is disabled, and with simulation NOFB, 
in which the emission of LW radiation is also disabled. 
In the last simulation, gas forms stars, but these stars 
do not emit radiation, and in addition, the intensity of 
the assumed LW background is set to zero. This sim- 
ul ation is there f ore id entical to simulation Z4 reported 
in iPawlik et al.l f|2011h , except for the inclusion of star 
formation. Important parameters of the simulations pre- 
sented here are summarized in Table [1] The final redshift 
of the simulations is z = 11. We note that simulations 
LW and NOFB are identical to the simulations used to 
es timate the detectabili ty of pair instability supernovae 
in lHummel et al.l (|2012f ). 

We use the friends-of-friends (FOF) halo finder, with 
linking paramet er b = 0.2, built into the substructure 
finder SUBFIND (|Springel et al.ll200fbl ). to extract halos 
from our simulations. Given a FOF halo, we use SUBFIND 
to identify its most bound particle and let it mark the 
halo center. We then obtain the virial radius, defined 
as the radius of the sphere centered on the most bound 
particle within which the average matter density is equal 
to 200 times the redshift-dependent critical density of the 
universe. The total mass inside this sphere defines the 
halo virial mass. We define the total SFR of a given halo 
as the sum of the SFRs of the gas particles it contains. 
Because we employ a stochastic star formation recipe, the 
rate at which star-forming gas particles are converted to 
star particles may randomly fluctuate around this SFR. 

We make use of the following relation between virial 
temperatur e T v -, r and virial mass M v -, r (e.g., Equ a- 
tion 3.12 in lLoebl [20101: see also iBarkana fc Loeb|[200ll) , 

4. FORMATION AND EVOLUTION OF THE DWARF 
GALAXY 

In this section we present the formation history of the 
simulated dwarf galaxy. We start with discussing the 
radiative feedback from LW and ionizing radiation on 
the assembly of gas and the formation of stars inside the 
dark matter halo hosting the galaxy (Section 14. ip . and 
on the assembly of the dark matter halo (Section I4.2|) . 
We finally proceed to investigate the robustness of the 



gaseous disks forming at the center of the dwarf galaxy 
halo under the radiative feedback (Section POl) . 

4.1. Baryon Assembly in the Dwarf Galaxy Halo 

Figure [5] shows the formation history of the dwarf 
galaxy in simulation LW+RT, including both photodis- 
sociating and photoionizing radiation (blue curves). For 
comparison, the figure also shows the formation histo- 
ries of the dwarf galaxy in the simulation in which star 
particles were sources of LW but not of photoionizing ra- 
diation (LW; red curves), and in the simulation in which 
star particles remained dark and the intensity of the LW 
background was set to zero (NOFB; black curves). The 
formation history is obtained by using SUBFIND to locate 
the dwarf halo progenitor that contains most of the 50 
most-bound particles of the dwarf halo at the final simu- 
lation redshift z = 11, and then repeating this procedure 
to find the progenitor of this progenitor and so on, trac- 
ing the halo assembly back to z — 25. The quantities 
displayed in Figure [5] are obtained from the properties of 
the particles inside the virial radius. The figure shows 
that the dark matter halo hosting the emerging dwarf 
galaxy grows by about three orders of magnitude in mass 
in about 300 millio n years, consistent with expectations 
(|Pawlik et alj|2011t their Figure 1). 

The evolution of the galaxy in simulation L W+ RT pro- 
ceeds in several main phases which will be discussed be- 
low: (1) the assembly of a dark matter minihalo with 
mass ~ 10 6 Mq and the accretion and condensation of 
gas inside it, leading to the formation of stars just below 
z = 23, (2) the subsequent accretion of gas under feed- 
back from star formation inside the dark matter mini- 
halo, (3) the evolution of this minihalo into an atom- 
ically cooling halo at z ~ 16, during which the prop- 
erties of the galaxy become insensitive to the inclusion 
of LW radiation, (4) the ensuing growth into a dwarf 
halo, during which the properties of the galaxy become 
robust against feedback from photoionization, and, fi- 
nally, (5) the formation of two nested rotationally sup- 
ported gaseous disks below z < 13.5. Our discussion of 
these phases will be complemented by comparisons with 
the evolution of the dwarf galaxy in simulations LW and 
NOFB. 

During the first phase, and as the dark matter halo 
grows in mass and accretes gas, both the baryon fraction, 
which is initially slightly smaller than the cosmic baryon 
fraction f2b/fi m ~ 0.17, and the central gas densities in- 
crease. By z = 25, the molecular hydrogen fraction has 
significantly departed from its initial value, enabling the 
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Fig. 1. — Gas density, temperature, ionized hydrogen fraction, and molecular hydrogen fraction (from left to right) averaged in a thin 
slice of thickness 0.2 kpc through the center of the minihalo progenitor at z = 22, shortly after the fir st stellar burst h as shut off. The 
relic photoionized region of diameter 4 — 6 kpc shows the typical non-spherical "butterfly" shape (e.g., Abel ct al. 1999). The dense gas 
inside the virial radius and along the fila ments cools rapidly, a llowing the residual free electrons to catalyze the formation of molecular 
hydrogen up to fractions r?H 2 ~ 2 X 10 — 3 (Oh & Haiman 2002). The molecular hydrogen fraction is also increased near the transition to 
the predominantly neutral IGM, a relict of the increased rate of molecular hydrogen formation inside the ionization fronts l lRicotti et alj 
2002). The dashed circle marks the virial radius. For purpose of visualization, this figure has been obtained from a simulation identical to 
simulation LW+RT bu t wit h an increased angular sampling, in which the propagation of photons was not limited to a single inter-particle 
distance. In Appendix IB. 31 we discuss that this change in parameters leaves the evolution of the minihalo essentially unaffected. 



minihalo gas to cool efficiently. The average gas temper- 
ature inside the virial radius is initially slightly higher 
than the virial temperature (dotted curve). Note that 
in simulation NOFB, as the molecular hydrogen fraction 
builds up, this relation then reverses and the virial tem- 
perature becomes higher than the average gas tempera- 
ture (e.g.. lO'Shea fe Normanll2007l ). At z « 23, the cen- 
tral gas densities become larger than the threshold den- 
sity for star formation. A single gas particle is turned 
into a star particle, triggering the emission of LW and 
ionizing radiation from the associated stellar burst. At 
that time, the halo has reached a mass of ~ 2 x 10 6 M . 

Photoionization from the radiation emitted by the first 
stellar burst almost instantly increases the average gas 
temperatures to > 10 4 K. The associated increase in 
thermal pressure pushes the gas away from the mini- 
halo center, and removes a fraction of it from inside the 
virial radius. As a consequence, the baryon fraction is 
reduced to about /bar = 0.15. The reduction in the 
baryon fraction is consistent with but slightly smaller 
than that found in previous simulations of the assembly 
of minihalos under feedback from st ar formation (e.g., 
I Wise fc Abeai2008at IWise et alJl2012T) . This may be be- 
cause the minihalo simulated here is fed by dense fila- 
ments, and the inflow of gas along these filaments pro- 
vides a strong obstacle for photoheating to drive the gas 
beyond the virial radius, and it reple nishes the photoe - 
vaporated regions with fresh gas (e.g. JAbel et al.ll2007D . 
Nevertheless, the reduction in the central gas mass is suf- 
ficient to drive the central gas densities below the thresh- 
old density for star formation, and the stellar burst shuts 
itself off. Figure [1] shows images of the gas density, tem- 
perature, and ionized and molecular hydrogen fraction 
around the minihalo at the end of the stellar burst. 

After the first stellar burst is shut off, the average 
mass-weighted molecular hydrogen fraction approaches 
~ 2 x 10 3 , as exp ected inside the relic H II region 
(Oh & Haiman 2002). The increased fraction of molec- 
ular hydrogen enables the gas to cool quickly, and the 
central gas densities increase. A few tens of Myr after 
the end of the first burst, the gas has become sufficiently 
cold and dense for ano t her stellar burst to be ignited 
(e.g.. lO'Shea et al.ll2005T: lAlvarez et al.ll2006l ). Photoion- 



ization heating from this second burst again lowers the 
gas densities, but this time there is no strong decrease in 
the baryon fraction. The accreting minihalo is thus mas- 
sive enough to retain most of the gas inside the virial re- 
gion. Because the negative feedback from photoheating 
is now less strong, this second starburst is more extended 
in time than the first one, and it involves the conversion 
of several gas particles to star particles. However, the 
combined feedback from the stellar clusters represented 
by the star particles eventually shuts off star formation. 
But already a few tens of Myr later the central gas densi- 
ties have again increased above the SF threshold density, 
and the galaxy continues to form stars. 

By redshift z ~ 16, the halo has reached virial tem- 
peratures T v j r > 10 K. Consequently, a significant frac- 
tion of the atomic hydrogen is collisionally excited, and 
its radiative de-excitation endows the gas with an ad- 
ditional channel to lose its thermal energy. Feedback 
from photoionization heating continues to keep the gas 
inside the halo at relatively low densities. However, the 
galaxy now forms stars continuously, albeit at a rate sig- 
nificantly smaller than in the absence of radiative feed- 
back. The comparison with the quickly rising SFRs in 
simulation LW that included emission of LW radiation 
but not that of ionizing photons shows that the feed- 
back from the photodissociation of molecular hydrogen 
by LW radiation alone becomes inefficient in preventing 
the gas from forming stars as the halo mass approaches 
the atomic cooling limit, in good agreement with pre- 
vious works (e.g.. lAhn fc Shapirol 120071: IWise fc Abel 
I2007al lO'Shea fc Norman! I2008T ). The transformation of 
the minihalo into an atomic cooling halo is accompanied 
by a significant decrease of the baryon fraction by « 20% 
at z ~ 15 — 16 and in the specific angular momentum of 
the gas j gas ee Jg as /Mg as , where J gas is the total gas 
angular momentum with respect to the motion of the 
most-bound particle and M gas the total gas mass. 

Below z < 15, photoionization heating becomes in- 
creasingly inefficient at reducing the density of the halo 
gas, and the gas then condenses up to > 10 5 cm -3 . As 
a result, the expansion of H II regions is impeded by 
the increased recombination rates. The properties of the 
galaxy then become rapidly insensitive to the inclusion 
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Fig. 2. — Assembly history of the dwarf galaxy. The blue solid curves in each panel show the evolution of the indicated property of the 
dwarf galaxy in simulation LW+RT which included both LW and ionizing radiation. For comparison, the red dash-dotted and the black 
dashed curves show, respectively, the evolution of the same properties in simulation LW that included only LW but not ionizing radiation 
and in simulation NOFB in which star particles remained dark and the intensity of the LW background was set to zero. The left-hand 
panel shows, from top to bottom: (a) virial mass M v i r (upper set of three curves which are nearly on top of each other) and stellar mass 
Mi, (lower set of three curves), and the mass corresponding to a virial temperature T v ; r = 10 4 K (dotted; Equation [6] with [i = 0.6); (b) 
SFR; (c) baryon fraction, and the universal baryon fraction Q^/flm ~ 0.17 (dotted); (d) specific angular momentum j ga s of the gas. The 
right-hand panel shows, from top to bottom: (a) the maximum gas density nu max! (b) tt> e mass-weighted mean temperature (T), and the 
virial temperature corresponding to the virial mass of the dwarf galaxy (in the NOFB simulation, dotted); (c) the mass- weighted mean 
ionized hydrogen fraction (rjHii), and the initial ionized hydrogen fraction, ?7hii = 3 x 10~ 4 (dotted); (d) the mass- weighted mean molecular 
hydrogen fraction (r)n 2 ), the initial molecular hydrogen fraction, 7Jh 2 = x 10 -6 (lower dotted), and the freeze-out molecular fraction 
expected in fossil H II regions, r)jj 2 ~ 2 X 10 — 3 (upper dotted; Oh & Haiman 2002). All quantities are computed considering only matter 
inside the redshift-dependent virial radius r v i r . The inclusion of photodissociation and photoheating implies as strong negative feedback 
on the condensation of gas and the formation of stars in the minihalo progenitor, but has little effect after the minihalo evolved into an 
atomically cooling galaxy. 



of ionizing radiation. The average mass-weighted ion- 
ized fraction, which has been increasing until then, re- 
mains roughly constant at around (rfan) 0.1. As a 
result of the large densities, the SFRs approach those 
seen in the simulation without radiation. The average 
mass-weighted molecular hydrogen fraction reaches val- 
ues in excess of (t7h 2 ) ^2 x 1CP 3 as H2 forms efficiently 
not only in fossil H II regions (e.g., IQh fc Haimai][200a 
Uohnson et al.l |2007|) but also in the ionization fronts 
of the H II regions unde r irradia tion from the ionizing 
stars (|Ricotti et al]|2002t see also. lShapiro & Kanj |l987[ 
Ueon et alll2012l ). 

A major merger at redshift z ~ 15 is accompanied by 
a strong increase in the specific angular momentum of 
the halo gas. The gas then settles in a disk, which is 
essentially in place at z ~ 13.5. The merger and the 
formation of this disk are shown in Figure 31 which dis- 



play a sequence of snapshots of the gas density inside 
the forming dwarf galaxy. The disk is surrounded by 
a second larger-scale disk at z w 12. The evolution 
and the properties of these two disks will be discussed 
in more detail in Section 14.31 below. After the forma- 
tion of the disks, the difference between virial and aver- 
age gas temperatures that has amplified since the halo 
has become an atomic cooling halo continues to increase. 
In the atomic cooling halo this difference arises primar- 
ily because a fraction of the gas does not shock-heat to 
the virial temperature upon accretion onto the halo. In- 
stead, the gas can efficiently cool inside the dense fila- 
ments that reach into the virial region an d supply the 
halo center with gas ()Pawlik et al.l |20TTI ). A qualita- 
tively similar difference between average and virial tem- 
peratures has been seen in ot her simulations of high- 
redshift low-mass galaxies (e.g-. IO'Shea fc Normanlr2007t 
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Fig. 3. — Dark matter density profiles of the dwarf galaxy in simulation LW+RT, which included both LW and ionizing radiation (bottom 
left; solid curves), and in simulation LW, which included only LW radiation (bottom right; solid curves). Curves of different colors show the 
density profile at different redshifts as indicated in the legend. The density profiles have been normalized by a singular isothermal density 
profile to reduce the dynamic range. In each panel, the set of dotted vertical lines on the left marks the redshift-dependent gravitational 
softening length e, and the set of dotted vertical lines on the right marks the redshift-dependent virial radius. For comparison, both bottom 
panels also show the dark matter density profile of the dwarf galaxy in the simulation without radiation (NOFB; dashed curves). The top 
panel in the left-hand (right-hand) figure shows the ratio of the dark matter density profile in the LW+RT simulation (LW simulation) and 
the dark matter density profile in simulation NOFB. The horizontal dotted line marks a ratio of 1, corresponding to the outcome in which 
the radiative feedback has no effect on the dark matter density profile. The inclusion of LW or ionizing radiation leads to a reduction of the 
central dark matter densities in the minihalo progenitor. As the minihalo turns into an atomically cooling galaxy and radiative feedback on 
the distribution of baryons becomes inefficient, the dark matter density profile approaches a singular isothermal profile in all simulations. 
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4.2. Effect of Radiative Feedback on the Dark Matter 

Halo 

Radiative feedback also affects the properties of the 
dark matter halo hosting the emerging dwarf galaxy. Fig- 
ure [3] compares the evolution of the dark matter density 
profile of the galaxy in simulation NOFB with that in 
simulation LW+RT (left panel), and with that in simula- 
tion LW (right panel). The bottom panels show the dark 
matter density profiles for each pair of simulations at four 
representative redshifts. The profiles are scaled by divid- 
ing by a singular isothermal profile pi so (r) — p vn r^ il: /r' 2 . 
In the last expression, r v ; r is the virial radius of the halo, 
p V i r = 200p cr it, and p cr it is the critical density. The top 
panels show the ratios of the density profiles of the sim- 
ulations including radiation and the simulation without 
radiation shown in the bottom panels, which helps to 
illustrate the effects of stellar feedback. In simulation 
NOFB, i.e., in the absence of radiative feedback (dashed 
curves), the dark matter density profile is approximately 
singular isothermal at all redshifts. 

The central dark matter densities in the simulations 
that included ionizing and/or LW radiation (solid curves 
in each of the bottom panels) are initially significantly 
lower, by up to factors ~ 5, than those in the simulation 
without radiation. The dark matter density profiles in 
these simulations thus do not follow a singular isothermal 
shape but show a spatially resolved dark matter "core" , 
extending to radii significantly larger than the gravita- 
tional softening scale (left-most vertical lines). The re- 



duction in the central dark matter densities is more dis- 
tinct and exists down to lower redshifts in simulation 
LW+RT than in simulation LW. The difference in the 
central dark matter densities originates in the difference 
in the distribution of the gas inside the assembling dwarf 
galaxy. In simulation LW, gas cannot cool and condense 
as efficiently as in simulation NOFB because molecu- 
lar hydrogen, the main coolant in low-mass primordial 
galaxies, is photodissociated by LW radiation. In simu- 
lation LW+RT, the central gas densities are, on average, 
further reduced as photoionization heating drives the gas 
away from the halo center. The radiative feedback on the 
distribution of baryons implies a significant change in the 
gravitational potential and, in turn, in the gravitational 
pull on the dark matter, which hence remains less cen- 
trally concentrated. 

The comparison of the dark matter density profiles 
in the simulations with and without feedback demon- 
strates that the ability of gas to cool and condense 
to high densities is crucial for establishing the singular 
isothermal den sity profile seen in the simulation withou t 
feedback (e.g., IWise fc Abel l2008bt iZemp et all [2012] ). 
That gas condensation can lead to a more concentrated 
dark matter density distribution is well known and is 
usually desc ribed using the framework of halo contrac- 
tion models dBlumenthal et al.l["l986HGnedin et al.ll200l 
iGnedin et all 120111 ). One may then expect the reverse 
of this process, i.e., the removal of gas from the halo 
center, to lead to a reduction in the concentration of 
the dark matter. Indeed, previous works have demon- 
strated the ability of SN explosions to lower the cen- 
tral dark matter densities i n dwarf galaxies with masses 
> 10 9 Mm at z < 10 (e.g.. iNayarro et al l [1991 iRicottil 
1 20031: IMashchenko etJtLl 120061 : iMashchenko et all 120081 : 
iGovernato et alj 120121 ). Here we have shown that ra- 
diative feedback can have a qualitatively similar effect 
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Fig. 4. — Assembly history of the disks in simulation LW+RT which includes both LW and ionizing radiation. The panels present face-on 
views of the assembling dwarf galaxy and show the evolution of gas densities in the redshift range 11.0 < z < 14.5 in cubical slices of 
linear size 0.5 kpc. The locations of star parti cles emitting LW an d ionizing radiation is indicated with black dots. This set of panels can 
be directly compared with that in Figure 8 of L Pawlik ct al. (2011), which shows the gas densities in a simulation identical to simulation 
LW+RT presented here except that it did not include star formation or radiation. Radiative feedback creates locally confined low-density 
cavities (best visible in the bottom right panel), but it is not sufficiently strong to prevent the assembly of the disks inside the atomically 
cooling dwarf galaxy. 



on the dark matter distribution in minihalos (see also, 
IWise fc AbeJl2008bh . 

4.3. The Disks 

Figures S] shows a sequence of snapshots of the gas den- 
sity distribution centered on the dwarf galaxy in simula- 
tion LW+RT. At the final simulation redshift, the gas at 
the halo center is organized in a massive central, spheri- 
cal clump, an inner compact disk, and an outer extended 
disk. The panels can be compared with Figure 8 of 
IPawlik et al.l (|201 ID , which shows the gas densities in our 
earlier simulation identical to simulation LW+RT dis- 
cussed here except that it did not include star formation 
or radiation. The comparison shows that the radiative 
feedback is weak and the inclusion of star formation and 
LW and ionizing radiation does not prevent the forma- 
tion of the disks seen in that earlier simulation. 

The sequence of events leading to the formation of 
the two disks is very similar with and without the in- 
cl usion of radiat i on, a nd the reader may therefore refer 
to IPawlik et al.l f|201 lh for additional details. A major 
merger at redshift z < 15 channels gas in the halo center. 
This leads to the formation of the first gaseous disk by 
z = 13.5. The disk subsequently develops spiral arms. 
The spiral arms exert torques that imply an outward 
transport of angular momentum in the disk gas. As a 
consequence, the disk shrinks in size. At z ~ 12, a se- 
quence of minor mergers replenishes the halo center with 
gas, producing the second gaseous disk. This disk re- 
mains spatially extended until the end of the simulation. 
The second disk surrounds the first disk, and the two 
have orientations tilted with respect to each other. This 
tilt is an interesting consequence of the hierarchical as- 



sembly of the emerging dwarf galaxy (e.g., iRoskar et all 
120101: (Pawlik et al.ll201ll) . The formation redshifts of the 
disks are similar in all simulations and hence insensitive 
to the inclusion of radiation. 

Photoheating creates low-density regions in the disk 
gas, leaving behind a disk morphology more complex 
than in the absence of photoheating (compare, e.g., the 
bottom right panel of Figure [H with th e corresponding 
panel in Figure 8 of IPawlik et aLl l20TlT) . These regions 
remain, however, locally confined, and the disks remain, 
as a whole, intact and relatively unaffected by photo- 
heating. This insensitivity of the gas to photoheating is 
consistent with the fact that at the time of formation of 
the disks, i.e., at z < 15, the galaxy halo has already en- 
tered the atomic cooling regime. Hence, its gravitational 
potential is sufficiently deep to confine photoheated gas 
in its interior. A positive feedback loop which consists 
of high gas densities confining the photoheated H II re- 
gions, allows the gas to collapse to increasingly higher 
densities, further confining the photoheated gas. 

The left panel of Figure [5] shows the Toomre param- 
eter, averaged in annular bins, of the disk gas at the 
final simulation redshift, Q = c s k/(7tGS), where c s is 
the adiabatic sound speed, k = (4S1 2 + rdfl 2 /dr) 1 / 2 the 
epicyclic frequency, an d f2 the angular velocity of the 
disk gas (|Toomrel ll964). In all three simulations Q > 1 
at radii larger than the gravitational softening length, 
implying that the disk configuration is stable against 
fragmentation. Note that the inclusion of LW and ion- 
izing radiation increases disk stability by increasing the 
sound speed. The inclusion of ionizing radiation further 
helps to preserve the disks because photoheating evapo- 
rates the gas from the low mass halos merging with the 
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Fig. 5. — Azimuthally averaged properties of the disks at z = 11 in the three simulations without radiation (NOFB; black dashed curves), 
LW radiation (LW; red dash-dotted curves), and LW and ionizing radiation (LW+RT; blue solid curves). Left: Toomre Q parameter of 
the disk gas. A value Q > 1 implies stability against gravitational fragmentation. Middle: SFR surface densities Ssfr as a function of 
distance r from the most bound particle. Right: Star formation rate surface densities Ssfr as a function of gas surface densities £ ga s- The 
dotted line shows the relation expected from the star formation recipe, Esfr = T^Sgas- The vertical dotted lines in the left and middle 
panels mark the softening scale e. 



dwarf galaxy. In the absence of feedback, on the other 
hand, baryon-rich low mass halos that pass through the 
disks can disturb t he disks significantly (e.g., Figure 3 in 
iPawlik et afll20Tl . 

The middle and right panels of Figure [5] show spheri- 
cally averaged profiles of the SFR surface densities (mid- 
dle), and the Kennicutt-Schmidt relation (right), i.e., the 
relation between SFR surface density and gas surface 
density, at the final simulation redshift. The Kennicutt- 
Schmidt relation is characterized by a power-law behav- 
ior and a suppression of star formation below gas surface 
densities < 10 2 M Q pc -2 . The suppression is a result 
of the star formation recipe, which limits star forma- 
tion to densities above 500 cm -3 . At the final simu- 
lation redshift, such densities are realized both in the 
central region and in the disks, and the galaxy shows 
a spatially extended morphology of star formation. In 
contrast, star formation in the galaxy before disk for- 
mation is limited to the central region (see Figure |4|. 
The Kennicutt-Schmidt power law behavior can be un- 
derstood by writing S S fr = Tg^Egas, where r gas = pjp± 
is the gas consumption time. Using r gas = (Equa- 
tion [I]), we find, Ssfr = t* Ems. This relation is shown 
with the dotted line in Figure [5j and the results from the 
simulation are in close agreement with it. 

Taken at face value, the large Toomre Q values that 
indicate stability against disk fragmentation seem incon- 
sistent with the nonzero star formation rates of the disks. 
However, our simulations do not have sufficient resolu- 
tion or the physical detail required to study fragmenta- 
tion of the disks into individual stars. We have therefore 
employed a phenomenological model for star formation 
according to which stars form from gas with densities 
larger than the adopted threshold density for star for- 
mation, an approach employed in most galaxy forma- 
tion simulations. It remains open if, at higher resolu- 
tion or greater physical detail, the disks in our simula- 
tion would fragment to form stars, or if they would re- 
main stable, possibly feeding a central massive black hole 
(e.g., lEisenstein fe LoeblH995t iKoushiappas et al J 120041 : 



lLodato fc Natarajanll2006f ). Addressing these issues in 
cosmological simulations such as the simulations here is 
computationally challenging. Simulations of isolated disk 
galaxies have shown that star formation is slower in disks 
that are more stable as quantified b y the smallest Toomre 
Q value in the disk (jLi et alJl2005t ). It may therefore be 
that star formation in high-redshift low-mass disk galax- 
ies is less efficient than implied by our simulations. 

5. REIONIZATION AND RADIATIVE FEEDBACK FROM 
THE FIRST STARS 

Figure |5] shows projections of the gas density (top) and 
temperature (bottom) in cubical slices through the re- 
finement region in simulation LW+RT at three represen- 
tative redshifts z — 19, 16, and 13 (from left to right). 
For comparison, the figure also shows the gas density 
and temperature at z = 13 in the corresponding slice 
through simulation NOFB (right- most panels). In the 
following, we discuss the impact of LW and ionizing ra- 
diation on the properties of the IGM and the formation of 
low-mass galaxies in the neighborhood of the simulated 
dwarf galaxy. 

The feedback processes exerted by radiation are well- 
known and are briefly summarized here. LW radiation 
photodissociates molecular hydrogen, reducing the abil- 
ity of low-mass halos to condense their gas. This sup- 
press es star formation, providing a negative feedback 
(e.g.. lHaiman et aT1ll997[ ). Ionizing radiation photoheats 
the IGM to ~ 10 4 K, and the implied Jeans filtering 
impedes the accretio n of gas into halos with virial tem- 
perat ures < 10 4 K ([Shapiro et alJll994t iGnedin fe Hul 
1998). In addition, photoheating evaporates gas from 
low-mass halos an d reduces the ab i lity of gas to cool 
insid e them (e.g.. lEfstathioul Il992t iThoul fc Weinberg! 
119961: iWiersma et all 12009ft . On the other hand, pho- 
toionization generates free electrons, which catalyze the 
formation of molecular hydrogen, thus increasing the 
abilit y of gas inside low-mass halos to cool an d form stars 
(e.g- lRicotti et al.l2002t \Oh fc rlaimanl2002ft . Photoion- 
ization therefore provides both a negative and a positive 
feedback on star formation. A comprehensive overview 
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Fig. 6. — Evolution of gas densities (top) and temperatures (bottom) in simulation LW+RT at three representative redshifts z = 19, 16, 
and 13 (from left to right). For comparison, the right-most panels show the gas densities and temperatures in simulation NOFB at z = 13. 
The cubical slices of fixed physical dimensions are centered on the comoving position of the most-bound particle in the dwarf galaxy at 
z = 10. The dashed circles are centered on the dwarf galaxy progenitor and have a radius equal to the virial radius of the progenitor. 
Photoheating raises the gas temperature to ~ 10 4 K. The implied increase in gas pressure evaporates gas from inside low-mass halos and 
smooths gas density fluctuations in the IGM. 



of the effects of radiativ e feedback can be found in, e.g., 
ICiardi fc Ferraral (l2005h . 

5.1. The Start of Reionization 

Figure [7J quantifies the ionization and thermal his- 
tory (left) as well as the evolution of LW intensities 
and of the molecular hydrogen fraction (right) around 
the dwarf galaxy in simulation LW+RT (blue curves). 
The neighborhood considered here is defined, for sim- 
plicity of geometry and to reduce boundary artifacts, as 
the sphere of comoving radius 5r v i r com (z = 10), where 
r V ir,com{z — 10) w 34 kpc comoving is the virial radius 
of the dwarf galaxy at z — 10. We center the sphere 
on the comoving position of the most-bound particle of 
this galaxy at z = 10, and this is the same position on 
which the slices in Figure |6] are centered. The spheri- 
cal neighborhood is contained in the refinement region 
at all redshifts. Our qualitative discussion below is not 
sensitive to the precise definition of the neighborhood 
adopted here. Volume-weighted (mass-weighted) aver- 
ages are obtained by weighting with the the SPH kernel 
volume (particle mass). 

The fraction of the volume ionized (solid) increases 
with decreasing redshift until z < 15, after which it re- 
mains approximately constant at (rjHii}v ~ 0.5. The 
fraction of the mass ionized (dashed) shows a qualita- 
tively similar behavior, but it slightly decreases after 
z < 15 and reaches (^Hii)m ~ 0.2 at the final simu- 
lation redshift. The redshift below which the ionized 
fractions no longer increase is similar to the redshift at 
which the galaxy evolves into an atomically cooling ob- 
ject. The fraction of the mass ionized is initially close to 
but slightly larger than the fraction of the volume ion- 
ized, but this changes at z < 19 after which the fraction 



of ionized mass becomes increasingly lower than the frac- 
tion of the ionized volume. Reionization thus proceeds 
from the inside-out, from the dense halo gas into t he dif- 
fuse IGM, as expected (e.g., iWise fc Abel l2008aft . The 
mass- and volume-weighted average gas temperature of 
the ionized gas, which we have defined here as gas with 
ionized hydrogen fraction ?7hii > 10" 3 , fluctuates be- 
tween 3000 and 10000 K, in good agreement with the 
results in Figures 12 and 13 in I Wise fc Abel (|2008al ). 

Figure [7] shows that the ionization of the region around 
the main galaxy in simulation LW+RT is accompanied 
by an increase in the average LW intensities. Before the 
formation of the first star, the average LW intensity is 
close to but slightly smaller than the intensity of the 
imposed LW background (dotted line), a consequence 
of self-shielding. Emission of stellar radiation increases 
the LW intensity above the background, and the latter 
quickly becomes unimportant. At z < 16 the LW inten- 
sities increase more rapidly and approach the intensities 
found in simulation LWthat only included LW radiation 
(red curves), and the LW radiation efficiently destroys 
the molecular hydrogen. However, the mass-weighted 
fraction remains significantly larger than the volume- 
weighted fraction, mostly because of self-shielding. The 
comparison with simulation LW shows that the inclusion 
of ionizing radia tion promotes the fo r mation of molecular 
hydrogen (e.g.. lHaiman et all 119961: iRicotti et all 120011: 
lOh fc Haimanll2002ft . 

The dwarf galaxy highly but not fully ionizes its neigh- 
borhood by z = 11. This is likely due to a combination of 
reasons. First, as the galaxy evolves into an atomically 
cooling object, photoheating is no longer able to substan- 
tially reduce the gas densities inside it (see Section [5]) . 
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Fig. 7. — Reionization history of the immediate neighborhood of the dwarf galaxy, defined by the sphere of comoving radius fa 100 kpc 
centered on the comoving position of the dwarf galaxy at z = 10, the same position used to center the images in Figure [6] In all panels, 
solid curves show volume- weighted averages, and dashed curves show mass-weighted averages. Top left: Mean ionized fractions (»?hii)- 
Bottom left: Mean gas temperatures (T). Thick (thin) lines show average temperatures including only gas that is ionized (neutral), i.e., 
gas with »?hii > 10~ 3 (rjjjn < 10 -3 ). Top right: Mean normalized LW intensities (J2i)- The average was computed in log-space to reduce 
the dynamic range, i.e. (J21) = 10^ log io ^21) . The dotted line shows the intensity of the imposed LW background. Bottom right: Mean 
molecular hydrogen fractions (??h 2 )- The dotted horizontal line marks the initial molecular hydrogen fraction, r)jj 2 = 1-1 X 10 — 6 . The dwarf 
galaxy highly but not fully reionizes the gas in its overdense neighborhood by the final simulation redshift, z = 11. 



Ionizing photons emitted by the stellar sources are then 
efficiently consumed by recombinations in the dense gas, 
thus limiting the ionizing impact on the surrounding 
IGM. Second, as we will discuss in Section lix^l below. star 
formation in the dwarf galaxy and its neighboring galax- 
ies exerts a strong negative feedback on neighboring low- 
mass galaxies, suppressing star formation inside them. 
These galaxies thus cannot significantly contribute to 
reionizing the gas. Finally, the considered volume lacks 
more massive galaxies, which would be robust against 
the feedback from the dwarf galaxy and could potentially 
help reionizing the gas. This is an artifact of the refine- 
ment region being too small to contain these rarer halos. 
Note also that the galaxy resides in a highly overdense 
region, and our computation of the mean ionized fraction 
includes contributions from the self-shielded neutral gas 
in halos inside this region. 

5.2. Feedback on Galaxy Formation 

Figure [8] shows the evolution of the minimum virial 
mass of star-forming halos in the refinement region, as 
found in simulation LW+RT hi which stars emitted both 
LW and ionizing radiation (LW+RT; blue). The cor- 
responding evolutions in the simulation in which stars 
emitted LW but not ionizing radiation (LW; red), and in 
the simulation in which no radiation was present (NOFB; 
black) , are also shown. We computed the minimum mass 
both by considering halos that form stars for the first 
time (crosses) , which is hereafter referred to as the mini- 
mum collapse mass, and by considering all halos, includ- 
ing those that have previously formed stars (boxes with 



matching but lighter colors). Our definition of the min- 
imum collapse mass is insensitive to our choice of the 
threshold density for star formation tih,sf = 500 cm~ 3 , 
at least in the absence of feedback from star formation, 
because the transition from densities nn > 10 cm~ 3 to 
n a ^ 10 3 cm" 3 which bracket the adopted star forma- 
tion threshold density then occurs within a narrow range 
of halo masses (see the bottom left panel in Figure ^j. 

In the absence of LW or ionizing radiation, the min- 
imum collapse mass is consistent with the mass of ha- 
los at virial temperature T vu - w 2200 K (Equation [5] 
with /i = 1.22) independent of redshift. This mass 
is consistent with but slightly larger than the masses 
~ 5 x 10 5 — 10 6 M Q at the time of gas collapse in pre- 
vious simulations of minihalos (e.g.. lYoshida et al.ll2003l: 
l O'Shea &: Normanl2007t iWise fc Abelll2008at iGreif et al l 
I2008T ) . The small difference is probably a numerical 
artifact of our limited resolution, which is lower than 
those realized in the previous works. Differences are also 
expected because the mass of the halo that forms the 
first star is subject to statistical uncertaintie s related to 
the l o w number of investigated halos (e.g., [Greif et al.l 
[200l lO'Shea fc Normanl l2008h . and because of differ- 
ences in the strength of dynamic al heating from tidal 
interactions and mass accretion (jYoshida et al.l 120031 : 
lO'Shea fc NoTnTanl 12001 . The inclusion of LW radi- 
ation increases the minimum collapse mass. This in- 
crease in the minimum collapse mass due to photodis- 
sociation by LW radiation has been in vestigated in de- 
tail in a number of previous w orks (e.g.. |Machacek et al.l 
l200lt IMesinger et all [200l lO'Shea fe Normanl 120081: 
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Fig. 8. — Minimum collapse mass, i.e., the mass of the lowest 
mass halo inside the refinement region that forms stars for the 
first time (crosses). The minimum collapse mass is shown both 
as inferred from simulation LW+RT (blue), which included both 
dissociating and ionizing radiation, and, for comparison, in the 
simulations without radiation (NOFB; black), and dissociating ra- 
diation only (LW; red). We also show the virial mass of the lowest 
mass halo that forms stars, including halos that have previously 
hosted star formation, in the same three simulations (squares with 
matching but lighter colors) . The dashed curve shows a virial mass 
corresponding to a virial temperature T vl[ = 2200 K (Equation [6] 
with /i = 1.22), which provides a good description of the evolution 
of the collapse mass in simulation NOFB. The minimum collapse 
mass is raised in simulation LW due to photodissociation of molec- 
ular hydrogen by LW photons and in simulation RT+LW addition- 
ally due to the Jeans filtering of the IGM. In the latter simulation, 
the enhanced molecular hydrogen fraction in fossil H II regions can 
lead to a reduction in the minimum collapse mass below that in 
simulation LW. 



ISafranek-Shrader et al]l20H . 

The inclusion of both LW and ionizing radiation in- 
creases the minimum collapse mass in a similar manner 
but often by a smaller factor, demonstrating a positive 
feedback from the enhanced formation of molecular hy- 
drogen in ionization fronts and fossil H II regions. Even- 
tually, however, the negative feedback from photohcating 
outweighs the positive feedback, and there is no new star- 
forming halo inside the refinement region below z < 13. 
Star formation can still proceed in halos that have pre- 
viously formed stars, but the masses of these halos are 
significantly larger than those in simulation NOFB. This 
is the result of the Jeans filtering of the IGM, which 
impedes the accretion of gas on low-mass halos. The 
increased scatter in the minimum collapse mass in the 
presence of photoionization is an expression of the local 
nature of feedback from photohcating, which is limited 
to inside the H II regions. 

The top panels of Figure [9] show the specific SFRs, 
i.e., the SFRs divided by virial mass, both as a func- 
tion of virial mass (left) and of stellar mass (middle), 
and the stellar mass fractions (right) for the halos in- 
side the high-resolution region at the final simulation 
rcdshift, z = 11. The bottom panels of Figure [9] show 
the maximum gas densities inside the virial radius (left) 



and the baryon mass fractions as function of virial (mid- 
dle) and stellar mass (right) for these halos. In the 
simulation without radiation (NOFB), the specific SFR 
~2x 10~ 10 yr" 1 is nearly independent of halo mass in 
the range *~ 10 7 - 10 9 Mq. The specific SFR is reduced 
in a fraction of the halos with masses < 10 7 M Q just 
above the minimum collapse mass, which we attribute 
primarily to the effects of dynamical heat ing during the 
gravi tational collapse of these halos (e.g., lYoshida et al.1 
2003). Indeed, the fact that the transition to high central 
gas densities occurs within a finite range of halo masses 
shows that dynamical heating plays a non-negligible role, 
and in the lowest mass halos this may be amplified by 
the limited mass resolution we afford. Note that some of 
the halos show an increased baryon fraction /b ar ^ 0.2 
as a result of dynamical interaction and ongoing mergers 
with other halos. 

The inclusion of LW radiation implies a complete sup- 
pression of star formation only in halos with masses 
< 2 x 10 7 Mq, corresponding to virial temperatures sig- 
nificantly below 10 4 K (vertical dashed line). This is 
in qualitative agreement with the results from previous 
high-resolution simulations of the collapse of minihalos in 
the presence of a LW radiation background. These works 
demonstrated that as the minihalo mass approaches the 
atomic cooling limit, the presence of LW radiation can- 
not prevent the build-up of molecular hydrogen, because 
it is catalyzed by the elevated electron f raction inside 
central structure formation sh ocks (e.g., iWise fc Abell 
I2007al lO'Shea fc Normanll2008[). and also becau se of self- 
shielding (e.g.. lAhn fc Shanirdl2007llSiisall2007l K In con- 
trast, the additional inclusion of ionizing radiation sup- 
presses star formation in essentially all halos below the 
atomic cooling limit by evaporating the gas from the halo 
centers. 

The simultaneous inclusion of both LW and ionizing 
radiation reduces the average baryon fractions in nearly 
the full range of simulated halo masses. The reduction 
of the baryon fraction is strongest, on average, for the 
lowest-mass halos with masses < 10 6 M , and for ha- 
los in the intermediate mass regime, with masses in the 
range 10 7 - 10 8 M . Because halo masses < 10 6 Mq 
are below the minimum collapse mass, halos in this mass 
range do not form stars, and hence their baryon frac- 
tion is reduced with respect to that in simulation NOFB 
due to the effects of radiation from external sources, and 
due to Jeans filtering in the photoheated IGM. More 
massive halos may accrete gas despite external feed- 
back and Jeans filtering, which explains why the baryon 
fraction increases, on average, in the halo mass range 
- 10 6 -10 7 M Q . Halos in the range 10 7 -10 8 M are effi- 
ciently forming stars, and therefore their baryon fractions 
are strongly red uced by radiative feedb ack from inter- 
nal sources (e.g.. iRicotti &: Gned in 2005). Finally, halos 
with masses > 10 8 M Q , corresponding to virial tempera- 
tures > 10 4 K, are robust not only against Jeans filtering 
and feedback from external sources, but also against the 
feedback from internal sources. The baryon fraction of 
these halos is reduced primarily because there was not 
yet enough time for accretion to compensate for the mass 
loss in past episodes of gaseous outflows. 

6. PROSPECTS FOR OBSERVATIONS WITH JWST 
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Fig. 9. — Properties of halos at z = 11, the final simulation redshift, in the refinement region in the three simulations without radiation 
(NOFB; black), with only LW radiation (LW; red), and with LW and ionizing radiation (LW+RT; blue). Each symbol represents a single 
FOF halo. Dashed and dotted vertical lines, if shown, mark the virial mass of a halo with virial temperature T v ] r = 10 4 K (Equation [6] 
with fi = 0.6), and the mass 100 jtidm- A horizontal dashed line, if shown, marks the universal baryon fraction. The top axis, if shown, 
labels the virial velocity v c = (GM vir /r^) 1 / 2 = 17.0 (M/10 8 M ) 1 / 3 [(l + z)/10] 1 ^ 2 km s" 1 (e.g., Equation 3.11 in lLoeb|[20Tol) . Top left: 
star formation rate, normalized by virial mass, as a function of virial mass. Middle left: SFR, normalized by virial mass, as a function 
of stellar mass. Top right: stellar mass fraction. The diagonal dotted line marks the stellar mass fraction implied by the presence of a 
single star particle. Bottom left: maximum gas density. The horizontal line marks the star formation threshold density. Bottom middle: 
baryon fraction as a function of virial mass. The curves of matching but lighter colors show medians to help guide the eye. Bottom right: 
baryon fraction as a function of stellar mass. All quantities have been computed considering only matter inside the virial radius r v ; r . 
Photodissociation of molecular hydrogen by LW radiation suppresses star formation completely only in minihalos with masses significantly 
below the atomic cooling limit. Photoheating, on the other hand, suppresses star formation also in more massive minihalos. The baryon 
fraction in simulation LW+RT is lowest in non-starforming halos with with masses < 10 6 Mq, in which it is strongly reduced by radiative 
feedback from external sources, and in star-forming halos with masses 10 7 — 10 s Mq, in which it is additionally reduced by radiative 
feedback from internal sources. 



In this section we estimate the fluxes of stellar and 
recombination radiation expected from the dwarf galaxy 
simulated here. We use these estimates to discuss the de- 
tecta bility of the first gala xies with the upcoming JWST 
(e.g.. lGardner et aT1l2006l) . We focus on the non- ionizing 
UV continuum at wavelengths around 1500 A (here- 
after UV1500), and on the Lya recombination line. We 
have present e d init ial estimates of the expected flux in 
iPawlik et all (|20T1 . The discussion here improves on 
our earlier work by basing the flux estimates on the dwarf 
galaxy simulation LW+RT presented in the current work. 
This simulation, which tracked the emission of LW and 
ionizing radiation from massive metal-free stars, provides 
us with direct predictions of the SFRs in low-mass high- 
redshift galaxies evolving under the radiative feedback 
from the first stars. 



JWST will image high-redshift galaxies in the UV1500 
continuum using NIRCam. JWST will also perform 
spectroscopic observations of these galaxies in the Lya 
line using NIRSpec. Observations of the UV con- 
tinuum are routinely used to infer t he SFR den- 
sity i n the high-redshift universe ('e.g.. iBunker et all 
2001; ISawicki fc Thompson] [2001 IBouwens et al.l 120091: 
Finkelstein et al.ll201 0I). a key quantity which not only 
constrains the nature of the stella r pop ulations (e.g., 
iDunlop et all l2012t IBouwens et~aH 120101) . but also al- 
lows one to address the capability of galaxies to reion- 
ize the universe (e.g..lStiavelli et al. 2004; IBouwens et al.l 
120091; iFinkelstein et al.ll2010D . Observations of the Lya 
line, on the other hand, have enabled, e.g., spectroscopic 
confirmation of gal axy candidates well into the epoch 
of reionization (e.g., iRhoads et all 120041; live et al.l 12001; 
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Fig. 10. — Flux of the combined stellar and nebular UV1500 continuum (left), and in the Lya recombination line (right) as implied by 
the SFRs of the dwarf galaxy in simulation LW+RT including LW and ionizing radiation (solid curves in the top panels, w hich addit i onally 
show the stellar mass in that galaxy with dotted curves and labels on the right-hand axis). The estimate is based on the Schacrcr (2003) 
stellar population synthesis models with zero metallicity and top-heavy IMF (red solid). For comparison, we also show the expected flux 
assuming metal- free stars with a normal IMF (dashed blue), and stars with metallicities Z = 5 X 10 — 4 Zq and normal IMF (dot-dashed 
black). The flux estimates assume that all the ionizing photons are absorbed inside the galaxy, i.e., / csc = 0. The line flux estimates can 
be rescaled to non-zero escape fractions by multiplication with (1 — / OS c). All flux estimates scale linearly with the SFR. Dotted horizontal 
lines show the sensitivity limits for observations with JWST, assuming exposures of 10 4 , 10 5 , and 10 6 s (top to bottom) and S/N=10 
UPanag ia 2005; Gardner ct al. 2006). Dwarf galaxies with SFRs ~ 0.1 Mq yr — 1 such as simulated here will be among the faintest galaxies 
JWST will detect in deep exposures of the z > 10 universe. 



iLehnert et aII[20Tot lOno et al.|[2012h . 

JWST will further enable spectroscopic observations 
of high-redshift galaxies in the Hel640 line, using NIR- 
Spec, and in the Ha line, using MIRI. The intrinsic 
strength of the Ha line is weake r by a factor of about 10 
than that of the Lya line (e.g., lSchaerer l l2003l ). Unlike 
the Lya line, the Ha lin e is, however, not affected by 
reson ant scattering (e.g., lLoeb fc Rvbickil I"l999t iSantosI 
j2004f ), rendering it a potentially competitive probe of 
high-redshift galaxy formation. The Hel640 line, on the 
other hand, is highly sensitive to the metallicity and 
IMF of the stellar populations, and a high ratio of lu- 
minosities in the Hel640 to Lya or Ha lines has been 
suggested a smoking gun for the existence o f massive 
meta l -free stars (e.g.. iTumlinson et al.l 120011 : lOh et all 
120011: iBromm et al.ll2001bO . A detection of the Hel640 
line would therefore put strong constraints on current 
theories of metal enrichment and star formation in the 
high-redshift universe, but is extremel y challenging be- 
cause of its weak intrinsic strength (e.g.. lZackrisson et al.l 
[Mil: iCai et al.|[201ll: IIn^l2lHll ). Because we find that 
both the Ha and Hcl640 line fluxes expected from dwarf 
galaxies such as simulated here are generally too weak to 
be observed with JWST, we do not further discuss these 
emission lines. However, the reader may refer to our ear- 
lier discussion of t he observability of th e first galaxies in 
Ha and Hel6 40 inlPawlik et all (|2011l) . 

We use the lSchaererf ()2003) population synthesis mod- 
els for constant star formation to convert the SFRs of 
the simulated dwarf galaxy in simulation LW+RT, which 
included both LW and ionizing radiation, into intrinsic 



luminosities in the L ya recombinati on line, and in the 
UV continuum. The ISchaererl ((2001 models require us 
to specify the IMF and the metallicities of the stellar 
populations, neither of which is predicted by our simula- 
tions and hence has to be assumed. To bracket plau- 
sible scenarios we repea t our analysis for three mod- 
els (see the discussion in IPawlik et al.l [20111 ). The first 
model assumes that stellar populations consists of metal- 
free ve ry massive stars , which we describe by employ- 
ing the lSchaere r (2003) model with zero metallicity and 
an IMF with Salpeter slope in the stellar mass range 
50 - 500 M Q (hereafter, top-heavy IMF). This is the 
same model employed to compute the ionizing and LW 
luminosities of the stellar bursts in our simulations. The 
second model assumes an IMF with Salpeter slope in 
the range 1 — 100 M (hereafter, normal IMF) and zero 
metallicity. The third and final model assumes the same 
IMF as the second model, but a non-zero (but low) metal- 
licity Z = 5 x 10~ 4 Z Q . 

We convert the line luminosities into observed line 
fluxes using 2 /(A ) = L(A e )/[47rd 2 j (z)], where L(X C ) is 
the intrinsic Lya line luminosity, and cZl(z) the luminos- 
ity distance to redshift z. The continuum luminosities 
are converted into observed flu x densities using an anal- 
ogous equation (Equation 5 in lPawlik et alll2011l) . The 
line luminosities and the nebular contribution to the UV 
continuum luminosities are derived assuming that all ion- 
izing photons are absorbed, i.e., that the escape fraction 

2 This differs from Pawlik ct al. (201]), where we converted the 
line flux into an equivalent flux density (Equation 4 in that work) . 
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is /esc = 0, thus maximizing the flux in the recombina- 
tion lines. The line luminosities can be rescaled to the 
case of non-zero escape fractions by multiplication with 
(1 — /esc)- We ignore the effects of resonant scattering 
on the observed Lya line fluxes. Lya RT simulations 
show that such effects can be very important, but they 
often depend sensitively on the specific structure of the 
invest igated galaxies and hence are difficult to generaliz e 
(e.g., iVerhamme et all 120081 ; iDiikstra fc Kramer! I2012T) . 
Finally, our flux estimates scale linearly with the SFRs. 

The curves in the bottom part of each panel in Fig- 
ure [10] show the resulting UV continuum fluxes (left), 
and the line flux in Lya (right) according to the three 
population synthesis models. For reference, the SFRs 
on which these flux estimates are based are displayed 
at the top of each panel (compare with the correspond- 
ing panel in Figure [2|) . We also show the total stellar 
masses at the top of each panel (right-hand axis). Fig- 
ure [TU] shows that the galaxy simulated here would be 
observable out to redshifts z < 13 in both the UV con- 
tinuum and the Lya line in deep surveys with exposures 
10 5 — 10 6 s. While the flux in the UV continuum is in- 
sensitive to the properties of the stellar populations, the 
Lya flux is significantly larger in the case of metal-free 
stellar population with top-heavy IMF than in the other 
two cases. At redshifts higher than z > 13, the reduction 
in the SFRs, mostly due to radiative feedback, causes the 
fluxes to decrease sharply below any practical detection 
limit. The resu lts presented here a re consistent with our 
earlier results (Pa wlik et al.l[20TTI ). after accounting for 
the differences in the SFRs. 

7. DISCUSSION 

The two nested gas disks in our simulations form only 
after the halo has reached a virial temperature signifi- 
cantly larger than 10 4 K. In light of the scale-free nature 
of cold dark matter structure formation, this relatively 
late formation of the disks may be surprising. However, 
it is consistent with results from previous zoomed simu- 
lations of the first atomically cooling galaxies which have 
not exhibited orderly ro tation of the gas (e.g. JWise et al.l 
120081 iGreif et al.ll2008D . On the other hand, the occur- 
rence of the disks in the dwarf-sized halos fits smoothly in 
line with simulations of galaxies more massive than the 
first atomic ally cooling galaxies assembling at lower red - 
shifts (e.g., IMashchenko et al.l 120081 iWise et al.l I2012D . 
Taken together, this suggests a threshold halo mass for 
the formation of the first disk galaxies in the range 
~ 10 8 - 10 10 M Q at z - 10. This mass scale is sim- 
ilar to the mass scale above which stellar feedback be- 
comes inefficient. However, the fact that simulations of 
the first atomically cooling galaxies have yielded a tur- 
bulent morphology even in the absence of star formation 
suggests that the two scales are physically unrelated. In- 
deed, in our simulations, the halo mass at the time of 
disk formation is insensitive to the inclusion of feedback. 
Nevertheless, the increased robustness of dwarf galax- 
ies ag ainst stellar feedback helps to preserve the disks 
(e.g., iKaufmann et al.ll2007f) . The prevalence of disks 
in halos with masses above 10 9 Mq that has recently 
been found in a simulation of galaxies in an overdense 
regio n at z ~ 10 probably mos tly owes to this latter ef- 
fect (|Romano-Di'az et al.ll20TTI ). 

Current observations of the classical dwarf satellite 



galaxies around the MW are consistent with a cored 
dark matter de nsity profile in at least a fraction of 
the o b jects (e.g..lGilmore et al.ll2007HJardel fc Gebhardl 
120121 iWolf fc Bullock! I2012D . These satellite galaxies 
are expected to have high-redshif t progenitors with to- 
tal m asses < 10 10 M Q at z > 6 (jBovlan-Kolchin et al.1 
|2012|) . similar to the total mass of the galaxy simulated 
here. Our results therefore support the hypothesis that 
such cored dark matter density profiles could originate 
in the gravitational interaction of the dark matter with 
the gas assembling in high-redshift low- mass halos un- 
der feedback from star for mation (e.g., iNavarro et al.l 
119961 iGovernato et al.| [2012) . Galaxies with a shallower 
dark matter density distribution may be expected to be 
more easily disrupted by tidal forces upon their entry 
in the MW halo and have their baryons stripped, thus 
possibly affecting the luminosities and observed abun- 
dance of the MW sate llites (e.g., IMashchenko &; Silisl 
l200l IMaverlrt~la1l2Tj07l) . 

Our results imply that galaxies with SFRs of ~ 
0.1 Mq yr _1 will be among the faintest galaxies JWST 
is like ly to detect at z > 10, confirming earlier estimates 
(e.g.. IHaiman fc Loeb! 119981 lOhl 119991: iTumlinson et al.l 
mOVL IZackrisson et al.l I20LU iPawlik et al.l 12011! ) . Ac- 
cording to our simulations, such galaxies reside in ha- 
los with masses of ~ 10 9 Mq having stellar masses of 
~ 10 7 M . In principle, JWST is sufficiently powerful to 
detect the light fro m stellar clusters with masses as low as 
10 5 - 10 6 MfT) (e.g I Johnson et~al~l l2009l: IZackrisson et al.1 
[20TTI: IPawlik et al l 12011ft . Our simulations do not sup- 
port the formation of such massive clusters, which would 
require local SFRs ~ 1 M Q yr _1 , an order of magnitude 
higher than found here. However, star formation in our 
simulations is unresolved, and it is possible that star for- 
mation in the first galaxies is more clustered or bursty 
than implied here. In this case, JWST could detect 
galaxies inside halos less massive than considered here, 
or inside halos that are more strongly affected by feed- 
back than suggested by our simulations. On the other 
hand, because we do not resolve the formation of stars 
from first principles, galaxies could be less efficient star- 
formers than implied by our simulations, and hence be 
fainter. Such fainter gala xies may still be seen if they are 
gravitational lensed (e.g., [Zackrisson ct al. 20T2h. 

Our simulations have ignored a potentially very impor- 
tant physical process, namely the explosion of massive 
stars in SNe. SNe can provide a strong negative feed- 
back by heating an d expelling gas from even relatively 
massive halos (e.g., iMac Low fc Ferraral 119991 ) . Previ- 
ous works have shown that SN feedback can suppress 
star formation strongly and lead to bursty star formation 
histories (e.g.. iStinson et aLll2007t ). SN feedback further 
may create a highly spatially inhomogeneous medium, 
likely enhancing the fraction of low column density sight- 
lines and hence the fracti on of esc aping ionizing ph otons 
(e.g.. lYaiima et al.l 120091 but see iDove et all 120000 . SN 
feedback may disturb the a ssembly of disks inside the 
first galaxies strongly (e.g.. IWise et al.l 120121 but see, 
e.g., IMashchenko et al.lT2008h . Moreover, SN feedback 
may modify th e structure of the dark matter halos sig- 
nificantly (e.g. .IMashchenko et al.ll2008l ; IGovernato et al.l 
l20TllBrooks fc Zolotovll2012D . 

Our simulations also ignored the chemical enrichment 
of the gas by the metals synthesized in stars. Sim- 
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ulations that track the production and transport of 
metals suggest that the transition between metal-free 
and metal-enriched stellar populations may occur early 
in the history of the universe (e.g.. iTprnatore et al.1 
[20071: IMaio et all l2010t IWise et alJ I2012h . Significant 
uncertainties, however, remain as to the efficiency of 
the mixing of metals with the primordial gas and the 
level of spatial homogen eity of metal enrichment (e.g., 
iScannapieco et al1l2002[ ). This, together with the fact 
that not all stars are expected to explode in SNe and en- 
rich the ga s but may instead collapse directly into black 
holes (e.g.. iHeger et al.l I2003T ) . leaves open the possibil- 
ity of the formation of metal-free stars in select regions 
of the universe down to relatively low redshifts (e.g. , 
" iTrenti et al.ll2009t I Johnson! [2010: 
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8. SUMMARY 

We have presented cosmological smoothed particle hy- 
drodynamics simulations of a dwarf galaxy assembling 
in a halo reaching 10 9 M at z — 10. The simulations 
were identical to our earlier simulations of such a galaxy 
in that they followed the non-equilibrium chemistry and 
cooling of primordial gas. They improved on our earlier 
simulations by including the formation of massive metal- 
free stars. To investigate the radiative feedback from 
these stars, we compared a simulation in which star par- 
ticles emitted both molecular hydrogen dissociating and 
hydrogen/helium ionizing radiation and a simulation in 
which star particles emitted only dissociating radiation 
with a simulation inside which star particles remained 
dark. 

Our main results are: 

• Dissociating and ionizing radiation exert a strong 
negative feedback by suppressing star formation in 
the main minihalo progenitor of the dwarf galaxy, 
but have little effect on star formation as soon as 
the progenitor evolves into an atomically cooling 
galaxy. 

• Radiative feedback suppresses the central dark 
matter densities in the dwarf galaxy main pro- 
genitor minihalo relative to the densities found in 
the simulation without radiation. The dark mat- 
ter density profile of the dwarf galaxy is singular 
isothermal independent of the inclusion of radia- 
tion shortly after the minihalo has evolved into an 
atomic cooling halo. 

• The dwarf galaxy halo hosts two nested disks below 



z < 12.5. The formation history and structure of 
the disks are insensitive to the inclusion of dissoci- 
ating and ionizing radiation. These results support 
a picture in which the first disk galaxies form in- 
side dark matter halos with masses > 10 9 Mq at 
z > 10. 



• The inclusion of dissociating and ionizing radiation 
lowers the baryon fractions inside the minihalos in 
the neighborhood of the dwarf galaxy. The baryon 
fractions are lowest in minihalos with masses < 
10 6 Mq, a consequence of Jeans filtering and pho- 
toevaporation from external ionizing sources, and 
in minihalos with masses ~ 10 — 10 s Mq, here pri- 
marily a consequence of photoevaporation of gas by 
internal ionizing sources. 

• Galaxies with star formation rates ~ 0.1 M Q yr _1 
will be among the faintest galaxies the upcoming 
James Webb Space Telescope will detect in deep 
exposures of the z > 10 universe. Our simu- 
lations suggest that such galaxies reside in halos 
with masses ~ 10 9 M Q and have stellar masses 
- 10 7 M . 

A major shortcoming of our simulations is the lack of 
feedback from supernova explosions. Such feedback can 
potentially have a significant impact on the evolution of 
low-mass galaxies. Feedback from supernovae may heav- 
ily disturb the assembly of disks, and strongly decrease 
the star formation rates inside dwarf galaxies. Our simu- 
lations also did not account for the chemical enrichment 
of the interstellar and intergalactic gas and the associated 
transition from metal-free to metal-enriched stars. The 
effects of supernova feedback and chemical enrichment 
are left to be investigated in future work. 
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APPENDIX 



In this appendix we discuss our implementation of photoionization and photoheating by stellar sources ( Appendix lA)) . 
We describe the im plem entation of the radiative transfer (RT; Appendix lA.lj) . the treatment of multi- frequency 
radiation (Appendix IA.2[) . the computation of the photoionization and photoheating rates (Appendix IA.3|) . and the 
coupling of the RT to the hydrodynamical evolution (Appendix IA.4[) . We also discuss tests of our RT implementation 
(Appendix iBl). specifically addressing the coupling of the RT with the solver for the chemical and thermal evolution 
(Appendix IB 1[) . and with the hydrodynamics (Appendix IB. 21) . We emphasize that the approximations and choices of 
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numerical parameters employed in this work and described below are specific to the current work. 

A. IMPLEMENTATION 

A.l. TRAPHIC 

The RT makes use of the RT code TR APHIC (IPawlik fc Schave 1 120081: IPawlik fc Schavd 120111) implemented in a 
customized copy of version 3 of gadget ()Springell I2005t ISchave et al.l I2010D . In simulations with traphic, the RT 
equation is solved by tracing a finite number of discrete photon packets emitted by ionizing source particles through 
the simulation box. This is done directly on the i rregular grid defined by the SPH particles, at the speed of light, and 
in a photon-conserving manner ( A bel et"aT][T99 9). In addition to the description of TRAPHIC given below, the reader 
may refer to the original publications for further details. 

Each photon packet carries photons of characteristic frequency v. We denote the total number of frequency bins 
used to discretize the radiation spectrum by N v . The transport of photons during a single time step proceeds by 
a succession of emission and transmission steps that move photon packets from individual particles, either SPH or 
star particles, to a number of N ng h SPH neighbors. The neighbors of a given particle are particles located inside the 
neighbor sphere, which is a sphere centered on the particle and contains N ng i, neighboring particles. Therefore, N ng \, 

determines the spatial resolution at which the RT is carried out. We choose N ng h = 48, which reflects a c ompromise 
between keeping high spatial resolution and controlling particle discreteness noise (|Pawlik fe Schave ll2008lh 

Star particles emit photon packets to their iV ng b neighboring SPH particles inside A^ec emission cones, each sub- 
tending a solid angle of 4-7T /N-^c ■ Emission cones tessellate the sky and are used to accomplish the isotropic emission 
of photon packet s to the SPH neighbors d espite the possibly highly anisotropic distribution of the SPH particles (see 
the Appendix in IPawlik fe Schavel 120081 ) . Photon packets emitted inside emission cones containing multiple SPH 
neighbors are split in proportion to the inverse squared distance between the source and the neighbors to account for 
the dilution of the radiation field with distance from the source. The central axes of the emission cones define the 
initial propagation directions of the emitted photon packets. The parameter A^ec determines the angular sampling of 
the sky as seen from any given ionizing source. We choose Nec = 8. Even though the orientation of the emission 
cone tessellation is randomly rotated at each RT time step, our choice of a small number of emission cones may imply 
increased random scatter in the distribution o f the photoionized and photoheated gas, especially around halos that 
contain only a few star particles (see Appendix IB. 31 for an illustration) . 

SPH particles that receive photon packets transmit these packets along the associated propagation directions to 
their downstream neighbors. Individual photon packets are transmitted only to the SPH neighbors located inside 
transmission cones, which are regular cones with solid angles subtending Att/Ntc steradians centered around the 
propagation directions and with their apex attached to the transmitting particle, where A^tc is a parameter specified 
below. The use of transmission cones prevents uncontrolled diffusion of photon packets on the set of irregularly 
distributed SPH particles and keeps the photon transport directed. Photon packets that are emitted inside transmission 
cones containing multiple SPH neighbors are split equally between these neighbors. 

The parameter A^tc sets the angular resolution at which the RT is performed. Because the transmission cones are 
defined locally at the positions of the transmitting SPH particles and because photon packets are only distributed 
among the subset of the N ng b neighbors that fall inside these transmission cones, the angular resolution is independent 
of the distance from the sources. In this work we adopt Nt c = 128. In test simulati ons of the RT around multiple 
sources inside a static cosmological density field presented in IPawlik fe Schavel (|2008l ) we found that the results had 
converged at a lower angular resolution of A^tc = 32. Note that virtual particles are created to accomplish the emissio n 
and transport of photon packets in cones that do not contain SPH neighbors, as described in IPawlik fc Schavel (|2011l ). 

Photon packets received by SPH particles are merged, which allows one to control the number of photon packets 
inside the simulation box and to avoid the scaling of the computational cost with the number of sources. The merging 
is done by binning photon packets in solid angle using a set of Nrc tessellating reception cones with solid angles 
47r/Afptc attached to each SPH particle. The merging is done separately for each frequency bin and hence it limits the 
number of photon packets to be transmitted next to at most Arc x N u . Accordingly, the computational cost of the 
RT scales with A^rc, but not with the n umber of ionizing sourc es. In this work we set N-rc = 8, motivated by test 
simulations similar to those presented in lPawlik~fc Schave (20QJ). 

Ph oton packets are associated a clock that is used to control the speed at which they travel (see IPawlik fc Schave 1 
2008). In this work we set this speed equal to the physical speed of light. However, for computational efficiency we 
allow photon packets to travel only a single inter-particle distance during individual RT time steps. In the limit of 
small RT time steps, transporting photon packets at the sp eed of light but limited b y a single inter-particle distance is 
equivalent to solving the time-dependent RT equation (see IPawlik fc Schave I [20081 for a similar argument in the case 
of the time- independent RT) . The hybrid approach employed here prevents the photons from traveling faster than the 
speed of light but, depending on the size of the RT time step and the inter-particle distance, may imply that photons 
travel at an effective speed that is lower than the speed of light, possibly resulting in an artificially delayed propagation 
of ionization fronts. We will discuss the effects of this approximation in Appendix IB.3I 

A. 2. Absorption of Ionizing Photons by Hydrogen /Helium 

A fraction 1 — exp[— t(z/)] of each photon packet in frequency bin v emitted or transmitted from particle i at 
position Vi to SPH neighbor j at position r^, separated by the propagation distance dij, is absorbed. The optical 
depth t{u) is the sum t(u) = ^r a {v) of the optical depths T a {v) of each absorbing species a e {HI, Hel, Hell}, 
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and r a {v) = f* 3 dr <7 a (v)n a (r) w a a (v)n a (rj)dij , The number of photons absorbed by species a is <5A" a b s , a (^) = 
KM/ S a Wai^jSNabs^), where the weights w Q (v) = T a (v) (IPawlik fc Schavd l20lTI) and 6N ahs (v) = £ Q SN^ hs , a (u) 
is the total number of absorbed photons. 

In this work we choose, for reasons of computational efficiency, to transport hydrogen-ionizing radiation using a 
single frequency bin, i.e., we set N v = 1. In this bin, we adopt frequency-averaged photoionization cross-sections (a^ n ) 
for the absorption by the individual species a using the grey approximation (e.g., Mihalas & Wcibcl Mihalas 1984), 



(ct 7Q ) = / dv — o ia [y) x 



dv- 



(Al) 



where J v (y) is the mean intensity of the radiation field, and hpv a is the ionization potential. We assume that the mean 
intensity J v (y) is characterized by a black body spectrum of temperature Tbb = 10 5 K, appropriate for the emission 
of radiation by the first stars (e.g. . lSchaererll2003l ) . which gives (ct 7 hi) = 1-63 x 10~ 18 cm 2 , (ct 7 h c i) = 2.28 x 10~ 18 cm 2 , 
and (q -YHeii) = 6.19 x 10~ 20 cm 2 . We have employed the fits to the photoionization cross-sections from lVerner et al.l 
(fl99l . 

At fixed mean intensity J v {v), the grey approximation implies photoionization rates identical to those inferred from 
a full multi-frequency treatment. However, because of the use of only a single frequency bin, our simulations ignore 
the hardening of the radiation spectrum with distance from the source caused b y the preferential absor ption of lower 
energy photons with larger absorption cross sections (see, e.g., the discussion in lPawlik fc Schavell201lD . 

A. 3. Photoionization and Photoheating Rates 

At the end of each RT time step the photoionization and photoheating rates are computed. We determine the 
photoionization rates directly from the total number of photons N a i, Sia (u) = ^A a b s ,a(^) absorbed by a given SPH 
particle during the RT time step At r in frequency bin v 1 



T 7Q = ^ , (A2) 



where Ah = "igas^H/niH is the number of hydrogen atoms associated with the SPH particle of mass m gas . This ensures 
that the same number of photons that have been removed from the simulation is use d to determine the ionization 
balance and temperature of the ionized gas, i.e., photon conservation (|Abel et al.lll999l ). 

The heating rate per atom due to photoionization of species a, assuming a single frequency bin in the grey approx- 
imation, is £ ia = T-y a (e a ), where 



(e a }= I dv ^ Jv ^ a ia (y)(h v v - h v v a ) x 



4nJ v (v) 
dv — (J ia (y) 

tlpV 



(A3) 



is the average energy of the absorbed ionizing photons in excess of the photoionization threshold. Because we assume 
that the mean intensity J v (y) is characterized by a black body spectrum of temperature Tbb = 10 5 K, the average 
excess energies for photoionization of hydrogen and helium are (em) = 6.32 eV, (e Hc i) = 8.70 eV, and (e Hc ii} = 7.88 eV. 

A. 4. Radiation- Hydrodynamical Coupling 

The radiation-hydrodynamical evolution of the gas is followed by invoking a series of subcycles to compute the 
dynamics, radiative transfer, chemistry, cooling, and heating of the gas. The dynamical evolution of the gas is 
followed on the gravito-hydrodynamical time steps set by the standard time integration scheme of the gadget code 
(Section l2.1[) . The RT is performed by subcycling the smallest (among all particles) gravito-hydrodynamical time step 
At s , with RT time steps of size At T . Unless stated otherwise, we adopt At r = min[Ai s ,0.1 Myr]. We assume that the 
chemical abundances and temperatures do not change during a single RT time step. 

Chemistry, heating, and cooling of the gas are computed at the end of each RT time step. This is done by subcycling 
the RT time step using time steps computed by the implicit non-equilibrium solver described in Section 12.21 Note 
that because photoionization rates are obtained from the number of absorptions computed under the assumption that 
species fractions and temperatures remain constant during the time Ai r , not all photons that have been absorbed 
always end up consumed in the computation of the chemistry and thermodynamics of the gas. This is a consequence 
of the evolution of the chemical composition and temperature during the subcycling. To ensure photon conservation we 
reinsert the remaining photons into the RT at the beginning of the next RT time step, and set them to be propagated 
along their original directions. 

The change in chemical abundances and temperature affects the dynamics of the particles only at the end of their 
gravito-hydrodynamical time steps, at which time the particle entropies are updated and the sizes of the new gravito- 
hydrodynamical time steps are determined. In gadget, particles evolve along a hierarchy of individual particle 
gravito-hydrodynamical time steps that may be much longer than the smallest such time step, the latter also being 
the time step at which the radiative cooling and heating is computed. This mean s that the radiation-hyd r odyna mical 
response of the gas to photoionization heating may occur with a delay. Recently, iDurier fc Dalla Vecchial (|2011l ) have 
shown, in the context of thermal feedback from supernova explosions, that the lack of a prompt response to localized 
energy injection may lead to a strong violation of the conservation of energy when using the standard gadget time 
integration scheme employed also here. Our simulations may suffer from similar numerical artifacts. 
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Fig. 11. — Test 4: cosmological rcionization by multiple stellar sources lllliev et al.l 120061) . Left: thin slice through the center of the 
simulation box at time t = 0.2 Myr, showing contours of neutral hydrogen fraction rfni = 0.5, on top of the density field (gr ey image). The 
ioniza tion fronts obtaine d in this work (red) are in very good agreemen t with the re ference results obtained using C2RAY (blue; Mcllcma ct al. 
2006) and CRASH (green; ICTardi et al.8 2001; Masclli ct al. 2003; Masclli ct al. 20Q3). Ri ght: histo grams of the g as temperature at t = 0.2 Myr. 
The results obtained in this work are in excellent agreement with those presented in Pawlik & Schayc (2011) (histogram labeled TRAPHIC 
thin for reference with our earlier work) obtained with a different solver for the chemical and thermal e volution. The y are in very goo d 
agreement with the reference solutions obtained with C2RAY, CRASH, and FTTE (Razoumov & Cardall 2005) as reported in Iliev et al, ( 2006), 
if one accounts for the differences in the atomic physics and the approximations employed. 



B. TESTS 

We have previo usly described a number of tests of our implementation of TRAPHIC in GADGET for the RT on static 
density fields. In iPawlik fc Schave I (f200l , we carried out monochromatic RT simulations of increasing complexity, 
assuming that the gas is compose d only of hydroge n and has a fixed temperature. By comparing with reference 
solutions such as those published in llliev et al.l ([20061 ). we demonstrated the ability of TRAPHIC to accurately capture 
the evolution of ionization fronts, and to reproduce the ionized fractions. We also showed that traphic is able to 
produce shar p shadows behin d op aque absorbers, which g ives rise to the typical "butterfly" shape of the ionized 
regions fe.g.. lAbel et~aT1 fl999h . In IPawlik fe Schavd (I2TjT1 . we extended our implementation of traphic to enable 
the transport of multi-frequency radiation, to account also for the ionization of helium, and to compute the evolution 
of the gas temperature due to radiative cooling and photoionization heating. Results fro m test simulations with this 
new implementation were in excellent agreement with reference solutions, such as those in llliev et al.l (12006D . 

The imp lementation of traphic in GADGET that we use in this work differs from that described in lPawlik fc Schaye I 
(2008) and IPawlik fc Schayc (201j]) in two main respects. First, we have substituted the explicit solver for the evolution 
of the chemistry and t emperature of primordial atomi c gas in the presen ce of photoionizations used and tested in 
IPawlik fc Schavd (|2011l ) with the implicit solver used in IPawlik et al.l (|2011[ ) and described in Section 12.21 This solver 
accounts for additional species, such as molecular hydrogen, and employs different rates for the atomic and molecular 
physics. It has been described and extensively tested in a number of publications (e.g., Uohnson" fc Bromml [20061 : 
iGreif et"aT1l2010l ). However, it has not yet been tested in combination w ith TRAPHIC . In Appen dix IB . 1 1 we therefore 
repeat the RT test simulation on a static cosmological density field from IPawlik fc Schavd ()201 lh . The results of this 
test are in excellent agreement with our previous results, demonstrating the validity of our implementation. 

Second, the current work is the first to employ traphic in radiation-hydrodynamic simulations that account for the 
feedback o f photo ionization heating on the gas dynamics. We have performed the radiation-hydrodynamical tests from 
llliev et al . (2009). and we have achieved excellent agreement with the published reference results in all of them. In 
Appendix IB . 2 1 we provide a discussion of one of the tests relevant to the present work, the simulation of the evaporation 
of a minihalo by an in terna l ionizing source. 

Finally, in Appendix IB. 31 we investigate how our adoption of a limited photon propagation speed and a finite angular 
sampling affects the early evolution of the minihalo in simulation LW+RT analyzed in the main text. 

B.l. Radiative Transfer and Chemical and Thermal Evolution 

In this section we re peat Test 4 of the cosm ological RT code comparison project (Ilie v et all 120061) that we have 
previously discussed in Pawlik fc Schave] ([201 ll ). using a chemistry solver different from the one employed here. The 
test involves the simulation of H II regions around multiple ionizing sources in a static cosmological density field. It 
was designed to capture important aspects of state-of-the-art simulations of hydrogen reionization, such as the delayed 




Fig. 12. — Test 6: profiles of the neutral and ionized fractions (left), temperature (middle), and gas dens ity (right) at t = 10 Myr. The 
results obtained with TRAPHIC (red solid) are compared with results obtained with FLASH-HC as published in Ilicv ct al. 2009 (blue dashed). 
The initial profiles at t = are also shown (black dotted). The agreement between the two results is very good. The small differences in 
the ionized fraction and temperature at r/0.8 kpc > 0.35 are primarily caused by differences in the numerical treatment of the blackbody 
spectrum. The monochromatic simulation with TRAPHIC employed a single frequency bin in the grey approximation, which cannot capture 
the pre-ionization and preheating found in the multi-frequency simulation with FLASH-HC. The shell of dense gas is slightly broader in 
the simulation with TRAPHIC than in the simulation with FLASH-HC, and its propagation is slightly delayed, leading to a slightly smaller 
radius at which the densities profile peaks. The small difference in the locations of the density profile peak explains the similarly small 
difference in the locations of the ionization front between the two simulations. The di fferences seen here are comparable to or smaller than 
the differences between the results obtained with a range of other RT codes shown in Iliev ct al. (2009). 



propagation in or trapping of ionization fronts by dense gas. 

The setup of this test is identical to that of Test 4 in IPawlik fc Schavel (|201 If ) . to which we refer the reader for a 
detailed description. Briefly, the initial conditions are provided by a snapshot (at redshift z rs 8.85) from a cosmological 
iV-body and gas-dynamical uniform-mesh simulation. The simulation box is Lbox = 0.5 hr 1 Mpc comoving on a side, 
where h = 0.7, and is uniformly divided into N cc n — 128 3 cells. We Monte Carlo sample this density field to replace 
the mesh cells with JVsph = N cc n = 128 3 SPH particles. The gas consists purely of atomic hydrogen and is assumed to 
be initially neutral at temperature T = 100 K. The ionizing sources are chosen to correspond to the 16 most massive 
halos in the box, and are assumed to have blac kbody spectra with temperature Tbb = 10 5 K. 

For comparison with|Pawlik & Schayc (201 lj) we solve the time-independent RT equation with an angular resolution 

of N c = 32, set the number of neighbors to which sources emit radiation to iV ng b = 32, employ a time step At r = 
10~ 4 Myr, and transport photons only over a single inter-particle distance per time step. We transport radiation using 
a single frequency bin, employing the grey photoionization cross-section (<7 7 hi) = 1-63 x 10~ 18 cm s^ 1 . We assume that 
each photoionization adds (em) = 6.32 cV to the thermal energy of the gas, as appropriate for the adopted blackbody 
spectrum. 

In Figure [TT] we show results of this simulatio n and compare th em with those presented in iPawlik fc Schavel (|2011l ) 
and also with the reference results reported in llliev et all (|2006[ ) . The agree ment is very good if on e accounts for 
the differences in the atomic physics and the approximations employed (see IPawlik fc Schavel 120111 for a detailed 
discussion). 

B.2. Radiation- Hydrodynamical Coupling 

In this section we perform Test 6 of the cosmological RT code comparison project (jlliev et all 120091 ) . This test 
consists of the simulation of the expansion of a H II region driven by an ionizing point source at the center of a 
spherically symmetric region with a steeply declining gas density profile. In contrast to the test described in the 
previous section, the gas distribution is not forced to remain static but may evolve in response to photoionization 
heating from the central source. The test situation shares some of the main features of the evolution of gas densities 
inside high-redshift cosmological minihalos under photoionization from a central massive metal-free star. 

The physical setup of the test is as follows. The initial gas densities are described by a cored isothermal density 
profile, 

n (r) = { U ° if r ~ r °' fBll 

n ^ r ' 1 no(r /ro)~ 2 otherwise. ^ ' 

We follow llTiev et all (|2009() and set no = 3.2 cm~ 3 and r$ = 91.5 pc, and we assume that the gas, consisting solely of 
atomic hydrogen, is initially neutral at temperature 100 K. The central source has a constant ionizing luminosity of 
10 50 phot ons s -1 throughou t the simulation, and is characterized by a blackbody spectrum with temperature 10 5 K. 
Following llliev etTa l (2009), all relevant cooling processes are included, except for Compton cooling. In this test, 
gravitational forces are ignored. 

We implement the power-law density profile by applying a local radial stretch along the direction towards the central 
source to the inter-particle distances of SPH particles that are initially distributed uniformly with density Denoting 
the initial particle positions with (f, 9, (f>), such a stretching is expressed by the transformation (f, 0) — > (r,9,<p), 
where 9 and <j> are the usual two angles needed to specify a position in spherical coordinates. Mass conservation then 
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Fig. 13. — Temperature averaged in a thin slice of thickness 0.2 kpc through the center of the minihalo progenitor at z = 22 in three 
simulations that are identical to simulation LW+RT discussed in the main text except for the treatment of the RT. Left: simulation 
test-base, with a single emission along JVec = 8 directions per RT time step and propagation at the speed of light but by at most 1 inter- 
particle distance, as in simulation LW+RT discussed in the main text. Middle: simulation test-speed, with a single emission along N^c = 8 
directions per RT time step and propagation at the speed of light (and no other limitation). Right: simulation test-sampling with 10 3 
emissions along JVec = 8 random directions per RT step and propagation at the speed of light (and no other limitation). The comparison 
of simulations test-base (left) and test-speed (middle) shows that the size of the final photohcatcd region is insensitive to the limit on the 
distance which photon packets can travel in a given RT time step, thanks to the photon-conserving nature of the RT with TRAPHIC. The 
comparison of simulations test-speed (middle) and test-sampling (right) shows that an increased angular sampling implemented by emitting 
photons along an increased number of random directions reduces artifacts in the shape of the photoheated region. Simulation test-sampling 
shown in the right panel has been used to generate Figure [l] 

requires the new coordinates (r, 0, <f>) to satisfy 

n H (r)r 2 sinQdrd6d4> = n^f 2 sin 6 dfdOdcj). (B2) 

Substituting a power-law density profile, n-a(r) oc r n , the last equation can be integrated to yield r oc f 3 /( 3 ~"). In the 
present case, n = 2, and hence r oc f 3 . Therefore, a singular isothermal profile can be obtained by stretching the initial 
inter-particle distances between uniformly distributed particles along the radial directions towards the center by a 
factor f 2 . The initial uniform density field is generated by placing particles randomly in the simulation box. To reduce 
the associated shot noise, the resulting particle distribution is regularized by evolving the particles under the influence 
of a reversed -sign (i.e., repulsive) gravitational force until the particles settle down into a glass-like quasi-equilibrium 
l|Whitel 1 1996ft . After generating the singular isothermal density profile, we replace the central region within ro with a 
uniform glass-like particle distribution to introduce the central core. 

The numerical parameters used in this test are as follows. We place the ionizing source in the center of a box with 
linear size 1.6 kpc. We sample the gas inside the core with radius using 30000 particles, which corresponds to an 
equivalent uni form grid resolu tion of 160 3 cells. This resolution is slightly higher than the uniform grid resolution 
employed in llliev et al.l (|2006ft . which was 128 3 cells. We use an angular resolution N c = 8, a number of neighbors 
for the emitting source N ng h = 32, employ a RT time step At T = 10~ 5 Myr, and transport photons only over a single 
inter-particle distance per time step. We limit the sizes of the individual particle hydrodynamical time steps to be 
less than 0.1 Myr. The radiation is transported using a single frequency bin and assuming a grey photoionization 
cross-section (<t 7 hi) = 1-63 x 10~ 18 cm s _1 , as well as that each photoionization adds (em) = 6.32 eV to the thermal 
energy of the gas. The parameter values are chosen to make contact with the conditions in the simulations presented 
in the main part of this paper, and to ensure numerical convergence. The test results are not critically sensitive to the 
adopted parameter values. However, we note that because of our choice of emitting photon packets only once per RT 
time step, the angular sampling depends on the RT time step, and a larger RT time step would imply an increased 
scatter in the ionized fraction, temperature and density profiles (see Appendix IB.3[) . 

The results of the test are shown and discussed in Figure [T2l The results obtained with TRAPHIC are in very good 
agreement with the reference results published in llliev et al.l (|2009l ) . 

B.3. Effect of Limited Photon Packet Propagation Speed and Finite Angular Sampling 

In this section we discuss the effects of our approximations regarding the limited photon propagation speed and the 
emission of photon packets by ionizing sources along a finite set of random directions per RT time step. To this end 
we have repeated the RT around the first stellar burst in the minihalo progenitor of simulation LW+RT. The first 
test simulation that we have performed, hereafter test-base, employed RT parameters identical to those employed in 
LW+RT. In particular, in this simulation ionizing sources emitted photon packets along iVsc = 8 random directions 
per RT time step, and the photon propagation was done at the speed of light but limited to a single inter-particle 
distance per RT time step. To investigate the effect of the limitation in the propagation speed, simulation test-base 
is compared with a simulation in which photon packets propagate at the speed of light with no additional limitation 
regarding the number of inter-particle distances per RT time step (hereafter test-speed; and with ionizing sources 
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emitting photons into Nec = 8 random directions per RT time step as in test-base). To investigate the effect of the 
angular sampling, simulation test-speed is compared with a simulation in which sources emit photon packets along 10 3 
times more random directions per RT time step (hereafter test-sampling; and transporting photons at the speed of 
light as in test-speed). The latter is implemented by repeating the emission of radiation along the Nec = 8 random 
directions 10 3 times per RT time step, with each emission using a different set of Nec = 8 random directions and 
emitting photon packets that contain by a factor 10 3 fewer photons. 

Figure [T51 shows the gas temperature in a thin slice through the center of the minihalo progenitor at z = 22, at the 
end of the stellar burst, in the three test simulations test-base (left), test-speed (middle), and test-sampling (right). 
The final extent of the photoheated region is similar in simulations test-base and test-speed despite the difference 
in propagation speeds. This is because the RT with traphic conserves the number of ionizing photons that are 
transmitted and absorbed, and because the final extent of photoheated H II regions is set primarily by this number. 
However, limitations of the speed at which photons are pr opagated can lead to dif ferences in extent and shape of the 
developing H II regions early during their evolution (e.g., iPawlik fc Schave Il2008|) . The increased angular sampling 
in simulation test-sampling reduces the noise that characterizes the shape of the photoheated region in simulation 
test-base, as expected. However, this noise is already small in simulation test-base, which hence yields a photoheated 
region similar in shape to that in simulation test-sampling. We note (but do not show) that the evolution of the 
properties of the minihalo, such as, e.g., the maximum gas density, is nearly identical in all three test simulations. 
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